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INTRODUCTION 


The  flow  field  in  the  aft-end  region  of  an  aerodynamic  vehicle  is  a  complex 
phenomenon  which  has  important  ramifications  for  the  vehicle  performance.  For 
example,  under  many  conditions  aft-body  drag  represents  a  significant  fraction 
of  the  total  drag  occurring  on  a  jet  engine  nacelle  and  this  portion  of  the  drag 
can  be  radically  influenced  by  the  nacelle  boattail  shape.  An  analysis  capable 
of  giving  accurate  predictions  of  the  flow  field  in  this  region,  including  off- 
design  operation,  could  play  an  important  role  in  the  boattail  design  process 
particularly  from  the  point  of  view  of  optimizing  the  trade-off  between  aft-end 
drag  and  boattail  weight.  The  aft-end  flow  field  problem  is  also  very  Important 
in  transonic  flow  situations  where  improved  design  can  lead  to  a  delayed  transonic 
drag  rise.  In  a  different  application  the  problem  also  becomes  significant  when 
considering  the  effect  of  a  propulsive  jet  plume  on  missile  aerodynamics.  In 
-this -ease -the_majx>r Jn.teres-t -centers,  .on  missile  stability  and  control  surface 
effectiveness.  Under  some  actual  flight  conditions  the  jet  plume  may  cause 
separation  and/or  significant  flow  field  changes  upstream  of  the  missile  base. 
These  in  turn  can  result  in  catastrophic  missile  pitch  up  or  loss  of  accuracy. 
However,  loss  of  control  surface  effectiveness  need  not  be  always  detrimental; 
it  can  be  an  advantage  in  the  transonic  portion  of  the  missile  boost  phase. 

Under  transonic  boost  phase  conditions  loss  of  effectiveness  may  lead  to  a 
decrease  in  missile  wind  sensitivity  which  is  a  desirable  effect  during  the 
boost  phase.  Thus,  in  the  case  of  a  missile,  loss  of  control  surface  effective¬ 
ness  during  the  transonic  boost  phase  combined  with  a  high  degree  of  control 
surface  effectiveness  in  supersonic  flow  may  be  the  desired  design  goal.  In 
addition,  the  recirculating  gases  can  of  course  pose  a  danger  because  of  their 
high  temperature.  Should  these  hot  gases  come  in  contact  with  an  unprotected 
surface,  catastrophic  failure  could  result. 

The  aft-end  flow  problem  has  been  the  subject  of  numerous  experimental  and 
theoretical  investigations.  Most  of  these  were  concerned  with  the  flow  field 
about  a  back  step  or  base  with  little  or  no  bleed.  However,  studies  more  germane 
to  the  present  problem  also  have  been  performed.  These  include  a  series  of 
experimental  efforts  aimed  at  investigating  the  missile  plume  effect  problem 
some  of  which  have  been  reported  in  the  bibliography  compiled  by  Batiuk  (Ref.  1), 
the  isolated  afterbody  problem  which  was  investigated  experimentally  by  Shrewsbury 
(Ref.  2)  and  the  jet  plume  problem  which  was  investigated  experimentally  by 
Strike  (Ref.  3)  and  Alplneri  and  Adams  (Ref.  4)  and  Galigher,  et  al  (Ref.  5). 
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Although  Refs.  1,  3  and  4  present  experiments  aimed  directly  at  the  jet 
plume  problem,  it  should  be  noted  that  accurate  wind  tunnel  simulations  of  the 
jet  plume  flow  field  are  difficult  to  attain  for  two  reasons.  The  first  problem 
which  arises  is  that  of  jet  exhaust  simulation.  In  general,  the  jet  exhaust  is 
represented  either  by  a  solid  body  jet  simulation  or  by  passing  either  hot  or 
cold  flow  through  a  model  exhaust.  At  present  it  is  an  open  question  as  to  how 
well  these  techniques  simulate  the  actual  flow  at  flight  Reynolds  numbers.  In  an 
investigation  of  missile  configurations  Burt  (Ref.  6)  concluded  that  his  plume 
simulators  give  approximately  the  same  afterbody  pressure  distributions  as  a  real 
plume.  Nevertheless  the  jet  simulation  problem  is  a  difficult  one  to  resolve  in 
general.  Secondly,  scale  effects  due  to  the  lack  of  simulation  of  both  the  aft- 
end  boundary  layer  and  the  free  mixing  layer  may  cause  uncertainty  in  the 
experimental  results.  This  is  particularly  true  in  the  transonic  regime  where 
viscous- inviscid  interacting  effects  play  a  major  role  in  determining  the  flow 
field,  and  where  the  interaction  effects  may  be  very  sensitive  to  the  viscous 
layer  thicknesses.  These  potential  problems  associated  with  wind  tunnel  testing 
show  that  theoretical  prediction  procedures  are  a  valuable  compliment  to  existing 
experimental  tools. 

As  in  the  case  of  experiments,  most  theoretical  analyses  of  the  aft-end  flow 
field  have  been  confined  to  back  steps  or  base  flows  without  bleed.  In  addition 
most  analyses  have  been  primarily  concerned  with  predicting  overall  flow  field 
quantities  such  as  base  pressure,  pressure  distribution  during  recompression, 
length  of  the  recirculation  region,  etc.,  rather  than  a  detailed  picture  of  the 
flow  field  behavior.  These  analyses  for  the  most  part  are  based  upon  'critical 
point'  methods  (e.g..  Refs.  7-11)  or  'component  analysis'  methods  (e.g..  Refs. 
12-15) .  Since  the  analyses  are  concerned  with  the  overall  rather  than  the 
detailed  nature  of  the  flow  field,  they  usually  have  relied  heavily  upon  boundary 
layer  Integral  methods,  mixing  analyses  and/or  semi-empirical  relationships. 
Nevertheless,  these  base  pressure  analyses  have  been  remarkably  successful  in 
predicting  certain  of  the  overall  flow  field  characteristics.  In  addition,  they 
are  usually  straightforward  to  implement  and  require  relatively  little  computer 
time.  In  supersonic  external  flow  where  the  major  concern  is  base  pressure  or 
centerline  recompression  pressure  rise  for  back  step  or  base  configurations 
with  simple  geometry  and  little  or  no  bleed,  these  rapid  techniques  may  be  quite 
adequate.  However,  the  inherent  quantities  which  make  the  procedures  rapid  and 
easy  to  implement  also  make  them  unsuitable  for  entirely  subsonic  flow  and  for 
predicting  details  of  recirculating  flow  and  make  them  questionable  for  cases 
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in  which  large  bleed,  complex  geometries,  etc.,  are  important.  For  example, 
these  methods  could  not  be  expected  to  predict  the  details  of  a  transonic  jet 
plume  interaction  flow  field  particularly  if  the  underexpanded  propulsive  jet 
were  to  cause  separation  of  the  boundary  layer  at  some  location  on  the  afterbody. 
Therefore,  if  a  detailed  description  of  the  aft-end  flow  is  desired,  an  alternative 
procedure  must  be  used. 

Recently,  the  advent  of  larger  and  faster  computers  and  more  efficient 
numerical  techniques  have  led  to  more  detailed  investigations  of  vehicle  aft-end 
flow  fields.  Several  investigations  have  used  the  Navier-Stokes  equations  to 
predict  backstep  and  base  flow  fields  (e.g..  Refs.  16-21),  and  only  recently 
have  the  Navier-Stokes  equations  have  not  generally  been  applied  to  the  full 
jet-plume  interaction  flow  field  problem  (Ref.  22).  In  addition,  the  subsonic 
boattail  problem  has  been  investigated  by  Rom  and  Bober  (Ref.  23)  who  used  an 
equivalent  body  to  represent  the  wake.  Fong  (Ref.  24)  combined  the  Chapman-Korst 
analysis  with  a  method  of  characteristics  analysis  to  study  the  jet  plume  problem 
which  included  the  effect  of  plume  induced  boundary  layer  separation.  Sinha, 

Zakkay  and  Erdos  (Ref.  25)  developed  an  invlscid  solution  to  the  underexpanded 
jet  plume.  More  recently  Grossman  and  Melnik  (Ref.  26)  modeled  the  transonic 
jet  plume  problem  by  interacting  three  separate  analyses.  The  jet  plume  is 
calculated  using  the  invlscid  method  of  Salas  (Ref.  27);  the  external  transonic 
flow  is  calculated  using  an  invlscid  procedure  based  upon  the  work  of  Murman 
and  Cole  (Ref.  28)  and  Jameson  (Ref.  29)  and  the  viscous  boundary  layer  and 
mixing  layer  are  calculated  using  the  integral  procedure  of  Green,  Weeks  and 
Broman  (Ref.  30).  Insofar  as  predicting  flow  details  and  interaction  effects, 
the  method  of  Ref.  26  represents  a  significant  improvement  over  earlier  'component 
analysis'  and  'critical  point*  methods.  (It  should  be  noted  that  these  latter 
methods  were  developed  primarily  to  predict  supersonic  base  pressure  and  recompres¬ 
sion  pressure  distributions  in  zero  and  small  bleed  cases  and  in  these  situations 
the  methods  seem  to  perform  quite  well.)  However,  as  Grossman  and  Melnik  noted 
in  Ref.  26,  their  analysis  does  make  compromises  in  the  treatment  of  viscous 
effects.  In  particular  the  method  uses  an  integral  boundary  layer  and  mixing 
length  analysis  and  such  integral  analyses  are  inherently  limited  by  their 
profile  assumptions.  Although  the  method  matches  viscous  and  invlscid  regions 
along  the  displacement  surface  in  regions  of  small  or  moderate  pressure  gradient, 
displacement  surface  matching  was  found  to  be  unstable  in  this  procedure  in 
regions  of  strong  pressure  gradient.  Therefore,  the  viscous  and  invlscid  flows 
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were  matched  along  curve  fits  joining  the  displacement  surfaces  on  either  side 
of  a  strong  pressure  gradient  region.  In  addition,  the  method  ignores  normal 
pressure  gradients  in  the  viscous  regions  and  does  not  rigorously  treat  either 
the  possible  boundary  layer  separation  on  the  boattail  or  regions  of  flow  reversal. 

In  the  present  report  we  focus  directly  on  the  region  of  the  flow  where  viscous- 
inviscid  interaction  are  important.  This  region  is  sketched  in  Fig.  1.  Consider 
the  interaction  of  an  under  expanded  propulsive  jet  with  a  supersonic  flow  over  a 
boattail  shaped  body.  The  Reynolds  number  is  assumed  sufficiently  high  so  that 
the  approach  boundary  layer  is  turbulent.  As  the  fluid  flows  over  the  boattail 
it  will  undergo  a  gradual  expansion,  up  to  some  point  upstream  of  the  boattail  base 
juncture  where  the  presence  of  the  jet  will  first  be  sensed.  The  location  of  this 
point  depends  in  general  upon  the  free  stream  and  jet  conditions  as  well  as  the 
geometry  of  the  boattail. 

In  the  region  of  upstream  influence  on  the  shoulder  of  the  boattail  the  flow 
will  compress  rapidly  forming  a  compression  wave  that  will  coalesce  into  a  shock. 

The  severe  adverse  pressure  gradient  that  is  generated  in  this  process  will  lead 
to  flow  separation  and  the  formation  a  recirculation  zone.  Immediately  downstream 
of  the  base,  the  flow  will  undergo  a  rapid  expansion,  and  thereafter  recompress 
with  the  formation  of  the  recompression  shock. 

Special  consideration  must  be  given  to  the  computational  domain  which  encompasses 
this  strong  interaction  region.  In  Fig.  2  a  schematic  of  the  computational  domain 
is  presented.  The  bounding  boundaries  of  the  region  of  interest  forms  an  L-shaped 
domain  containing  a  sharp  reentrant  corner  at  point  h. 

The  present  effort  considers  only  supersonic  inflow  and  supersonic  outflow 
boundaries  for  simplicity  and  since  the  amount  of  upstream  influence  in  the  (under¬ 
expanded)  propulsive  jet  is  expected  to  be  minor  the  upstream  conditions  for  the 
calculation  within  the  nozzle  are  applied  at  the  nozzle  exit  plane  (line  a-b  of 
Fig.  2) .  Consideration  of  subsonic  inflow  and  subsonic  exit  conditions  is 
conceptually  simple  but  probably  would  require  a  larger  domain  and  could  be  the 
subject  of  a  later  investigation.  The  upstream  conditions  for  the  flow  over  the 
boattail  are  applied  upstream  of  the  nozzle  exit  plane  so  as  to  allow  calculation 
of  possible  boattail  separation  effects  (line  c-d) .  The  outer  boundary  conditions 
for  the  external  flow  is  taken  along  a  line  parallel  to  the  nozzle  centerline 
which  is  far  enough  removed  from  the  nozzle  so  as  to  allow  application  of  freestream 
boundary  conditions  (line  d-e)  and  keep  spurious  shock  and/or  Mach  reflections 
from  the  boundary  from  impinging  on  regions  of  interest.  Finally,  the  outflow 
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boundary  of  the  present  region  of  interest  is  taken  downstream  of  the  recirculation 
zone  and  weak  streamwise  derivative  boundary  conditions  are  applied  along  this 
line  (line  e-f ) .  The  solution  at  the  downstream  boundary  could  be  used  as 
inflow  boundary  conditions  for  a  conventional  viscous  or  inviscid  procedure  which 
could  be  used  to  predict  the  downstream  "far"  flow  field. 

Since  the  computatonal  domain  is  L-shaped, standard  solution  algorithms  must 
be  modified  to  allow  for  efficient  computations.  Furthermore,  the  reentrant  sharp 
corner  is  a  geometrically  singularity  in  the  flow  field,  and  must,  therefore,  be 
treated  appropriately.  A  fuller  discussion  of  these  matters  will  be  given  later 
in  this  report. 

Under  the  present  effort  a  local  flow  field  analysis  is  considered.  This 
detailed  solution  then  could  be  joined  either  to  inviscid  solutions  or  viscous 
forward  marching  solutions,  as  appropriate,  in  an  interactive  manner  to  determine 
the  entire  flow  field.  For  supersonic  bounding  flows  the  matching  of  the  detailed 
Navier-Stokes  solution  to  the  external  flow  is  quite  straightforward,  provided 
of  course  that  all  regions  of  upstream  influence  of  downstream  disturbances  are 
contained  within  the  Navier-Stokes  solution  domain.  For  subsonic  or  transonic 
bounding  flows  the  transient  interacting  technique  of  Briley  and  McDonald  (Ref.  31) 
would  be  used,  either  as  a  linearized  correction  to  the  external  stream  or  as  a 
matching  surface  for  a  time  or  relaxation  step  in  an  inviscid  flow  calculation. 

The  procedure  that  is  used  in  this  effort  is  the  Multidimensional  Implicit 
Nonlinear  Time-Dependent  (MINT)  procedure  of  Briley  and  McDonald  (Ref.  32).  This 
represents  an  accurate,  efficient  numerical  procedure  for  solving  the  Navier- 
Stokes  equations.  The  procedure  was  originally  developed  for  laminar  duct  flow 
and  was  subsequently  applied  to  turbulent  duct  flow  by  Briley,  McDonald  and 
Gibeling  (Ref.  33).  It  has  since  been  applied  to  the  three-dimensional  combustor 
flow  problem  by  Gibeling,  McDonald  and  Briley  (Ref.  34).  This  procedure  represents 
a  well-proven  and  well-exercised  method  of  solving  the  Navier-Stokes  equations. 

The  MINT  procedure  has  been  previously  described  in  Refs.  32  and  33,  and  con¬ 
sequently  only  a  brief  outline  will  be  given  here  and  the  scheme  presented  in  detail 
in  Appendix  B.  The  method  can  be  briefly  outlined  as  follows:  The  governing 
equations  are  replaced  by  an  implicit  time  difference  approximation,  optionally 
a  backward  difference  or  Crank-Nicolson  scheme.  Terms  involving  nonlinearities 
at  the  implicit  time  level  are  linearized  by  Taylor  expansion  about  the  solution 
at  the  known  time  level,  and  spatial  difference  approximations  are  introduced. 

The  result  is  a  system  of  multidimensional  coupled  (but  linear)  difference  equa¬ 
tions  for  the  dependent  variables  at  the  unknown  or  implicit  time  level.  To 
solve  these  difference  equations,  the  Douglas-Gunn  (Ref.  35)  procedure  for 
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generating  alternating  direction  implicit  (ADI)  schemes  as  perturbations  of 
fundamental  implicit  difference  schemes  is  introduced.  This  technique  leads  to 
systems  of  coupled  linear  difference  equations  having  narrow  block-banded  matrix 
structures  which  can  be  solved  efficiently  by  standard  block-elimination  methods. 

The  method  centers  around  the  use  of  a  formal  linearization  technique 
adapted  for  the  integration  of  initial-value  problems.  The  linearization 
technique,  which  requires  an  implicit  solution  procedure,  permits  the  solution 
of  coupled  nonlinear  equations  in  one  space  dimension  (to  the  requisite  degree 
of  accuracy)  by  a  one-step  noniterative  scheme.  Since  no  iteration  is  required 
to  compute  the  solution  for  a  single  time  step,  and  since  only  moderate  effort  is 
required  for  solution  of  the  implicit  difference  equations,  the  method  is  computa¬ 
tionally  efficient;  this  efficiency  is  retained  for  multidimensional  problems  by 
using  ADI  techniques.  The  method  is  also  economical  in  terms  of  computer  storage, 
in  its  present  form  requiring  only  two  time-levels  of  storage  for  each  dependent 
variable.  Furthermore,  the  ADI  technique  reduces  multidimensional  problems  to 
sequences  of  calculations  which  are  one-dimensional  in  the  sense  that  easily-solved 
narrow  block-banded  matrices  associated  with  one-dimensional  rows  of  grid  points 
are  produced.  Consequently,  only  these  onc-dimensional  problems  require  rapid- 
acess  storage  at  any  given  stage  of  the  solution  procedure,  and  the  remaining  flow 
variables  can  be  saved  on  auxiliary  storage  devices  if  desired.  A  more  complete 
description  of  the  linearization  technique,  the  alternating  direction  technique 
and  the  solution  procedure  is  presented  in  Appendix  B. 
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THEORETICAL  ANALYSIS 


Governing  Equations 

In  this  section  the  governing  equations  are  presented  in  vector  form  for  a 
two-dimensional  flow.  The  resulting  Jacobian  transformed  governing  conservation 
equations  in  cylindrical  polar  coordinates  are  given  in  Appendix  A.  It  is 
convenient  for  turbulent  flows  to  formally  eliminate  density  fluctuations.  This 
can  be  accomplished  by  using  the  mass  weighted  or  Favre  average  (Ref.  36).  The 
governing  equations  are  then  identical  to  the  laminar  equations  with  the  flow 
variables  being  taken  as  mean  quantities  and  viscosity  being  taken  as  the  sum  of 
the  molecular  viscosity  p,  and  the  turbulent  viscosity  p^.  The  resulting  equations 
are 

Continuity 

(i) 


Momentum 


dp  u 

"aT 


+  V  -  (/3uu  )  =  -Vp  +  V  '  (  7 r  +  7TT) 


(2) 


=  T 

In  the  above  equation  n  and  n  are  the  average  and  turbulent  stress  tensor 
respectively. 


Energy  Equation 

In  the  present  formulation  it  is  desirable  to  write  the  energy  equation 
in  terms  of  the  static  enthalpy  h  because  of  simplifications  in  the  turbulence 
time  averaging,  viz. 


+  V  •  (pu  h)  =  -  V- (q  +  qT)  +  ^-  +  $  +  p*  ^ 

where  $  is  the  mean  flow  dissipation  defined  in  Eq.  (10)  and  e  is  the  turbulence 
kinetic  energy  dissipation  rate  (cf.  Eq.  (20)). 
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Turbulence  Kinetic  Energy  Equation 

In  the  present  work  a  turbulence  kinetic  energy  equation  and  a  specified 
length  scale  equation  have  been  utilized.  Following  the  derivations  of  (Refs. 
37  -  38)  and  Ref.  (39),  one  obtains 

<3(  /ok)  _  u. 

— £ —  +  V-{  ouk)  -  V-(  V  k )  W 

dt  cr  k 


+  /iT  [2D:D  -  -|-(7-u)2]  -  -|-/>kV-U  - 


and  k  is  the  turbulence  kinetic  energy 


k  -  ~  u'.u' 


Stress  Tensors 


The  stress  tensor  assuming  a  Newtonian  fluid  is 


=  2p.ID  -  (-§-/!  -  Kb)v-UII 


where  is  the  bulk  viscosity  coefficient  and  is  the  deformation  tensor,  i.e. 


[  ( vu)  +  (vO)T] 


The  turbulent  flow  stress  tensor  is  modeled  using  an  isotropic  eddy  viscosity 
formulation,  i.e., 


*  -p  u'u'  =  2^-tD  -  —  ( /iTV-  U  +  pk  )E 


where  k  is  the  turbulence  kinetic  energy,  Eq.  (5).  The  turbulent  viscosity 
must  be  determined  using  a  suitable  turbulence  model. 

Heat  Flux  Vectors 

=  -+t 

The  mean  heat  flux  vector  q  and  the  turbulent  heat  flux  vector  q  may  be 

written  as 


q  *  -  <[  Vt] 

qT  -  -*T  [  vf] 


(9) 

(10) 
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T 

where  ”  is  the  mean  thermal  conductiveity,  and  k  is  a  turbulent  conductivity, 
Mean  Flow  Dissipation 

The  mean  flow  dissipation  term  appearing  in  the  energy  equation,  Eq.  (3), 
is  defined  as 

*  *  2/aD  :d  Kb)(7-u)2  (11) 

In  the  present  analysis  Eqs.  (1-4)  are  solved  in  conjunction  with  Eqs.  (6-10) 

=  T 

and  the  constitutive  relations  for  y,  pT>  K^,  k  and  k  ,  which  are  given  in  the 
subsequent  section. 

Constitutive  Relations 

The  necessary  constitutive  relations  include  an  equation  of  state,  a 
caloric  equation  of  state,  a  turbulence  length  scale  distribution,  the  molecular 
viscosity  and  thermal  conductivity, 

Equation  of  State 

The  equation  of  state  is  that  of  a  perfect  gas 

p  =  /jRT  0-2> 


where  R  is  the  gas  constant, 

The  caloric  equation  of  state  is  taken  as 


e 


cvT 


(13) 
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where  would  be  dependent  on  the  gas  composition,  but  not  temperature.  We  may 
also  define  the  static  enthalpy  as 


h  .  CpT 


(14) 


where  c^,  the  specific  heat  at  constant  pressure  is  not  a  function  of  temperature. 
In  the  calculations  presented  in  this  report  Cp  and  c^  were  taken  as  constants. 

Molecular  Viscosity,  Bulk  Viscosity,  and  Thermal  Conductivity 

The  molecular  viscosity  for  the  gas  is  determined  from  Sutherland's  law. 


JL  .  I  I  lalh 

Mo  '  To  /  T+S, 


(15) 


where  =>  110°K  for  air. 

The  bulk  viscosity  will  be  assumed  to  be  zero  at  present, 


K 


B 


»  0 


(16) 


The  thermal  conductivity  may  be  determined  from  a  relation  similar  to 
Sutherland's  law,  e.g., 


2  T0+S2 

<0  'To'  T  +  S2 


where  S2  “  194 °K  for  air. 
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Turbulence  Model  Relations 

The  turbulent  viscosity  introduced  in  Eq.  (8)  is  obtained  from  the 
Prandtl-Kolmogorov  relation,  viz. 

2 


/*t"  C 


pk 

H-  € 


(18) 


and  the  dissipation  rate  is  given  by 


€  •  C, 


_  3/2 
3/4  k 


(19) 


where  the  turbulence  length  scale,  l,  must  be  specified  consistent  with  the 
expected  turbulence  structure  in  the  flow.  Following  Ref.  AO  the  constant  o^ 
will  be  taken  as 


°k  =  '-° 


(20) 


There  are  two  options  available  for  modeling  the  turbulence  near  a  wall. 

In  the  first,  grid  point  resolution  normal  to  the  wall  must  be  sufficient  to  define 
the  viscous  sublayer,  in  which  case  the  boundary  conditions  are  relatively  straight¬ 
forward.  However,  the  difficulty  with  this  approach  is  that  the  physics  of  low 
Reynolds  number  (transitional)  turbulence  must  be  modeled  in  a  reasonable  manner 
by  the  governing  turbulence  equations  (e.g.,  Jones  and  Launder,  Ref.  40).  An 
alternative  approach  is  to  employ  a  less  refined  mesh  and  force  the  turbulence 
variables  to  yield  values  consistent  with  a  boundary  layer  wall  function  formulation 
at  the  first  grid  point  away  from  the  wall.  The  difficulty  with  this  approach  is 
that  the  validity  of  the  wall  function  formulation  is  questionable  under  the  flow 
conditions  of  the  viscous  inviscid  interaction  problem.  A  transition  model,  which 
was  successfully  used  by  Shamroth  and  Gibeling  (Ref.  41)  in  a  time  dependent 
airfoil  flow  field  analysis,  has  been  implemented  in  the  computer  code.  The 
analysis  of  Ref.  41  follows  the  integral  turbulence  energy  procedure  of  Refs.  42  -  44 
by  utilizing  a  turbulence  function  a^,  where 


o,  *  CM,/2  /2 


(21) 


and  a^  is  taken  as  a  function  of  a  turbulence  Reynolds  number  of  the  form 
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a 


t 


1.0  +  6.66  a. 


t(RT) 

100 


(22) 


where 


a0  »  .0115 

f(Rr)  -  100.  RT°  a  RT  <  I 
fC  Rt )  -  68  I  Rt  +  614.3  RT>40 


(23) 


and  a  cubic  curve  was  fit  for  values  of  R^  between  1  and  40.  As  previously 
discussed.  Ref.  42  -  44  utilized  an  integral  form  of  the  turbulence  kinetic 
energy  and,  therefore,  R^  was  defined  as  an  average  value. 


Rr  ■  T-C'W  T  l  *  >'t,y* 


(24) 


In  the  present  effort  R^  was  defined  as  the  local  ratio  of  turbulent  to  laminar 
viscosity,  a-^  was  evaluated  via  Eq.  22  and  related  to  a^  via  Eq.  21. 

The  effective  viscosity  is  defined  as  the  sum  of  the  laminar  and  turbulent 
viscosities 


Meff  =  T 


(25) 


The  effective  conductivity  will  be  modeled  using  an  effective  Prandtl  number 
obtained  from  knowledge  of  turbulent  flows  of  gases  and  gas  mixtures,  i.e.. 


eff 


K  +  K 


cp/xeff 

Preff 


(26) 


and  Pr  *  0.9  for  air. 
eff 
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Turbulence  Length  Scale 

The  near  wall  region  mixing  length  is  obtained  from  the  McDonald's  model 
(Ref.  45)  with  Van  Driest  damping  (Ref.  46), 

-p- ]  [  i-exp(-y+/A+)]  <27> 

*-ao  J 

where  k  is  the  von  Karman  constant  and  A+  is  the  van  Driest  damping  coefficient, 


k  =  0.4 
A*  =  26.0 

and  1  *  .096. 

00 

The  nondimens ional  distance  y+  is  defined  as 


and  the  friction  velocity  u^  in  the  present  analysis  is  taken  as 


(28) 


UT 


(29) 


where  the  local  shear  stress  is  obtained  from 

ta  =  (  2D:  P),,z  (30) 

Note  that  for  small  y  the  tanh  function  in  Eq.  (27)  reduces  to  ky  while  for 
large  y  it  approaches 

A  more  detailed  description  of  the  actual  turbulence  models  used  for  the 
computation  and  their  effect  on  the  computation  will  be  given  in  the  discussion 
section  of  the  paper. 
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COORDINATE  TRANSFORMATION 

The  set  of  governing  partial  differential  equations  which  model  the  physical 
processes  was  presented  in  the  previous  section.  For  generality  these  equations 
were  written  in  vector  notation;  however,  before  these  equations  can  be  incorporated 
into  a  computer  code,  a  coordinate  system  must  be  chosen.  The  governing  equations 
can  then  be  cast  in  a  form  reflecting  the  choice  of  the  coordinate  system.  The 
coordinate  system  for  a  nozzle  boat-tail  geometry  must  allow  for  the  boat-tail 
body  to  be  a  coordinate  line,  the  nozzle  exit  plane  to  be  a  straight  line  and  it 
should  have  the  flexibility  of  clustering  points  in  regions  of  large  shear. 

Another  feature  of  the  coordinate  system  should  be  its  ability  to  adapt  the  mesh 
distribution  as  regions  of  large  shear  develop  during  the  calculation,  i.e. 
be  time  dependent.  In  order  to  permit  consistent  calculation  of  the  Jacobian  in 
a  moving  coordinate  system,  the  governing  equations  should  be  transformed  with  a 
Jacobian  transformation  of  the  form 

y’  *  yJ(x,,  x2,7j,i)  (31) 

T  =  t 

where  (x  ,  x2>  x.3)  are  the  original  coordinates  (cf.  Ref.  47).  In  cylindrical 
polar  coordinates  (x^  would  correspond  to  (r,  0,  z).  The  velocity  components 
remain  the  components,  (U^,  U2>  0^)  in  the  (x^,  x2>  x^)  coordinate  directions 
respectively.  The  new  independent  variables  yj  are  the  computational  coordinates 
in  the  transformed  system.  The  coordinate  system  requirements  for  the  problem  under 
consideration  may  be  represented  by  a  subset  of  the  general  transformation,  Eq.  (31) 


y1  =  y'u,,x3,t) 

yz  =  y2(x2)  (32) 

y3  =  y3(x,,  x3, t) 
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which  is  a  general  axisymmetric  time-dependent  transformation.  For  the  nozzle 

and  boattail  geometry  which  is  axisymmetric,  Eq.  (32)  reduces  to  y2  -  JL  and 

2  2 
all  derivatives  3/ 3y  are  assumed  to  be  zero. 

Application  of  the  Jacobian  transformation  requires  expansion  of  the  temporal 

and  spatial  derivatives  using  the  chain  rule,  i.e.. 


d<ft  deft  £  j  d<ft 
- +.?  y.'  dyl 


at 


dr 


)*• 


and 


dxi  j=i  *  dyJ 


(33) 


(34) 


where 


4  = 


dyJ 

dt 

dyj 

dx  i 


(35) 


The  relations  Eq.  (33-35)  are  first  substituted  into  the  governing  equations 
(1-4  )  written  in  Cartesian  or  cylindrical  polar  coordinates.  Then  the  resulting 
equations  are  multiplied  by  the  Jacobian  determinant  of  the  inverse  transformation. 


J  = 


d(x,  ,x2,  x3) 

d(y',y2,y3) 


dx. 

dx, 

dxi 

dy1 

dy2 

dy3 

dx2 

o# 

xl 

ro 

O' 

X 

ro 

dy  1 

dy2 

dy3 

dx3 

dx3 

d^3 

dy' 

dy2 

dy3 

(36) 


and  the  equations  are  cast  into  a  "semi-strong"  conservation  form  (Ref.  56) 
using  the  following  relations, 


3 

1 

i=* 


d(Jy,\)  _ 


dyJ 


=  0 


(37) 
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•VWM> 


ii  +  y  3jylt 

dr  j=i  dyj 


The  semi-strong  conservation  form  implies  that  all  factors  involving  the  radial 
coordinate  r  *  x^  remain  as  they  were  before  the  Jacobian  transformation.  The 
resulting  equations  are  presented  in  Appendix  B. 

The  geometric  relations  Eq.  (37-38)  may  be  obtained  from  the  transformation 
relations  for  Jy,^  and  Jy,|  in  terms  of  the  inverse  transformation  derivatives 
(e.g..  Ref.  44), 


Jy',l  =  *2,2  X3,3  *2,3  *3,2 
Jy',2  =  *3,2  *1,3  "  *3,3  *1,2 


J^,3  ’  X!,2  *  2,3  X  1,3  X  2,2 


Jy'  =  *2,3  ?3,  '  X  2,1  *3,3 


JyZo  =  x, ,  x,  i  -  x, -  *, , 
3  ,2  3,3  1,1  3,1  1,3 


Jy%  =  *1,3  *2,.  -  *.,.  *2,3 
Jy3,,  =  *2,.  *3,2  "  *2,2*3,. 
Jy3,2=  *3,.  V  -  *3,2*.,. 

Jy%  =  *.,.  *2,2  "  *1,2  *2,1 


Jyjf  =  -  £jyi 

w  y’k  aT 


Grid  Transformation 


The  flow  over  a  nozzle  boattail  configuration  entails  regions  of  large 
shear  near  the  surface  of  the  boattail  and  in  the  shear  layer  coming  off  the 
trailing  edge  as  seen  in  Fig.  2,  curve  h-g.  We  must,  therefore,  allow  for  mesh 
clustering  in  both  the  streamwise  and  radial  direction  at  the  body's  surface  and 
in  the  interior  of  the  flow  field.  In  addition,  the  boattail  surface  should 
correspond  to  a  coordinate  line  in  order  to  enable  one  to  easily  impose  the 
appropriate  surface  boundary  conditions. 

The  methods  employed  in  the  construction  of  the  coordinate  system  are 
based  on  generalizations  of  the  Roberts  type  transformation.  Ref.  (48).  Several 
different  analytic  forms  were  used  to  cluster  the  grid  points. 

In  the  computational  space  the  grid  spacing  is  uniform  and  for  simplicity 
set  to  1,  i.e.  Ax  *  1.0.  In  the  present  problem  the  constraints  of  the  grid 
transformation  are  chosen  so  that  the  boundaries  of  the  computational  domain 
fall  on  grid  points.  In  addition,  the  boundary  of  the  boattail  surface  is 
chosen  to  lie  on  grid  point  number  m,  for  instance.  The  choice  of  the  grid 
spacing  between  points  m  +  1  and  m  is  based  on  physical  considerations  such  as 
the  requirement  that  the  first  point  above  the  surface  should  lie  within  the 
sublayer  of  the  turbulent  boundary  layer.  An  analogous  procedure  is  used  to 
distribute  the  grid  in  the  streamwise  direction,  but  in  this  case  the  base  is 
fixed  at  a  desired  grid  point. 

This  type  of  procedure  gives  the  desired  degree  of  flexibility  for  the 
problems  under  consideration.  Furthermore,  since  the  grid  distribution  in  the 
streamwise  and  radial  directions  are  independent,  the  grid  will  in  general 
be  nonorthogonal.  This  feature,  rather  than  being  an  impedement  turns  out  to 
be  beneficial  in  that  we  do  not  overly  constrain  the  problem  to  conform  to  say 
an  orthogonality  condition,  which  would  necessarily  distort  the  boattail  base 
(cf.  Ref.  22)  when  using  a  body  fitted  coordinate  system.  By  allowing  for  a 
nonorthogonal  grid  the  shape  of  the  boattail  can  be  preserved.  Obviously,  the 
sharp  corner  at  the  Intersection  of  the  base  and  shoulder  introduces  a  singular 
point  into  the  calcultion  which  must  be  treated  appropriately.  In  the  following 
section,  this  matter  is  taken  up  and  discussed. 

Treatment  of  the  Corner  Points 

Corner  points  such  as  point  h  of  Fig.  2  where  the  boattail  shoulder  intersects 
the  base  present  special  problems  to  finite  difference  schemes  and  this  flow 
region  has  been  the  subject  of  several  recent  investigations.  An  overview  of  the 
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problem  is  given  by  Roache  Ref.  (49)  who  points  out  that  this  region  may  be  j 

subject  to  large  truncation  error.  In  addition,  as  discussed  by  Roache  the  1 

question  of  whether  the  potentially  large  truncation  error  is  confined  to  the  j 

immediate  region  of  the  corner  or  whether  it  has  a  much  larger  global  effect  is  ] 

still  unresolved.  ■ 

Several  types  of  methods  have  been  proposed  for  treating  this  region. 

A  common  practice  is  to  eliminate  the  corner  by  rounding  it,  and  thus  avoid  the 
problem  completely.  But,  this  method  changes  the  geometry  under  consideration. 
Alternatively  at  the  corner  a  local  analysis  can  be  undertaken  Refs.  (50-51). 

As  discussed  by  Fox  (Ref.  50),  Fox  and  Sankar  (Ref.  51)  represent  the  singular 
behavior  of  the  flow  by  using  special  equations  based  upon  a  series  expansion 
in  the  immediate  vicinity  of  a  wall  singularity  (such  as  the  wall  corner  point) . 

A  second  approach  which  also  isolates  this  area  has  been  developed  by  Ladeveze 
and  Peyret  (Ref.  52)  who  match  an  asymptotic  solution  of  the  flow  in  the  vicinity 
of  the  corner  with  a  finite  difference  solution  of  the  Navier-Stokes  equations 
elsewhere.  But  this  approach  leads  to  a  more  complicated  calculation,  since 
the  embedded  local  analysis  must  be  matched  to  the  main  flow  field  computation. 

The  inherent  errors  such  as  interpolating  onto  the  "large"  domain  could  offset  the  gains 
of  using  an  embedded  domain  analysis  technique.  Furthermore,  the  two  domain 
algorithm  assumes  that  the  grid  is  sufficiently  fine  in  the  corner  region.  This 
may  lead  to  undue  strain  on  the  grid  spacing  in  the  "large"  domain  or  additional 
grid  points  must  be  included,  leading  to  increased  run  times. 

As  a  compromise  between  these  two  approaches  the  "double  valued"  technique 
was  chosen.  As  discussed  by  Roache  (Ref.  49)  the  method  to  be  utilized  is  based 
upon  discontinuous  values  of  dependent  variables  at  the  corner  point.  As  shown 
in  Fig.  2,  the  corner  point,  h,  can  be  approached  from  either  a  horizontal  or 
vertical  direction.  If  the  required  boundary  condition  is  a  derivative  condition, 
the  value  at  the  corner  will  depend  upon  whether  the  derivative  is  taken  in  the 
horizontal  or  vertical  direction.  The  situation  is  resolved  by  taking  h  as  a 
double-valued  point;  i.e.,  a  second  grid  point  is  located  a  minute  distance  from 
the  corner  and  dependent  variables  may  be  discontinuous  between  the  two  points 
representing  the  corner.  One  of  the  points  would  be  used  in  evaluating  derivatives 
in  the  horizontal  direction  and  the  other  point  would  be  used  in  evaluating 
derivatives  in  the  vertical  direction. 

In  this  technique  the  corner  remains  sharp,  maintaining  the  original  geometry, 
but  the  singular  point  is  allowed  to  exist  in  that  the  values  of  the  pressure  and 
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density  are  permitted  to  vary  depending  on  the  coordinate  direction  from  which 
the  corner  is  approached.  This  type  of  approach  fits  in  naturally  with  the  ADI 
procedure  since  the  double  values  are  associated  with  the  appropriate  implicit 
direction.  It  is  important  to  note  that  only  quantities  whose  functional  values 
are  not  known  at  the  corner  assume  a  double  value,  i.e.  non  function  boundary 
conditions  are  applied  there.  Since  the  no-slip  condition  is  applied  at  the  wall, 
the  velocity  components  are  single  valued  and  are  set  equal  to  zero.  However, 
the  pressure  and  density  become  double  valued  at  the  corner  since  the  boundary 
condition  used  to  evaluate  these  variables  involve  derivatives.  If  a  gradient 
condition  were  used  for  the  wall  enthalpy  (heat  transfer  rate  prescribed)  then 
the  enthalpy  would  also  be  double  valued  at  the  corner.  In  our  case,  the  wall 
enthalpy  is  prescribed  so  that  a  single  value  of  the  enthalpy  was  employed. 
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SOLUTION  PROCEDURE 


The  computer  code  is  based  upon  an  axisymmetric  version  of  the  highly 
efficient  consistently  split,  linearized  block- implicit  solution  procedure 
(MINT)  for  the  compressible  Navier-Stokes  equations  developed  by  Briley  and 
McDonald  (Refs.  53,  32,  33),  and  subsequently  extended  to  multi-component, 
chemically  reacting,  turbulent  flows  by  Gibeling,  McDonald  and  Briley 
(Ref.  34).  This  procedure  solves  the  Navier-Stokes  equations  written  in  pri¬ 
mitive  variables;  in  the  MINT  procedure,  the  governing  equations  are  replaced 
by  either  a  Crank-Nicholson  or  a  backward  time  difference  approximation.  Terms 
involving  nonlinearities  at  the  implicit  time  level  are  linearized  by  Taylor 
series  expansion  about  the  known  time  level,  and  spatial  difference  approxima¬ 
tions  are  introduced.  The  result  is  a  system  of  two-dimensional  coupled 
linear  difference  equations  for  the  dependent  variables  at  tht  unknown  or 
implicit  time  level.  These  equations  are  solved  by  the  Douglas-Gunn  (Ref.  35) 
procedure  for  generating  ADI  schemes  as  perturbations  to  fundamental  implicit 
difference  schemes.  This  technique  leads  to  systems  of  one-dimensional  coupled 
linear  difference  equations  which  are  solved  by  standard  block-elimination 
methods,  with  no  iteration  required  to  compute  the  solution  for  a  given  time  step. 
An  artificial  dissipation  term  based  upon  either  a  cell  Reynolds  number  criterion 
or  the  rate  of  change  of  the  dependent  variable  may  be  introduced  selectively  into 
the  scheme  to  allow  calculations  to  be  performed  at  high  local  values  of  the  cell 
Reynolds  number. 

The  use  of  an  implicit  solution  procedure  requires  that  equation  coupling 
and  linarization  be  considered.  Both  of  these  topics  are  reviewed  in  detail 
by  McDonald  and  Briley  (Ref.  54)  and  Briley  and  McDonald  (Ref.  32).  There  it  is 
shown  that  for  a  given  spatial  grid  the  errors  arising  from  time  linearization 
of  the  nonlinear  terms  at  the  unknown  time  level  should  be  no  greater  than  the 
discretization  errors.  Weinberg  and  McDonald  in  Ref.  (55)  also  demonstrated 
this  result  for  a  model  problem.  Furthermore,  reduction  of  the  time  step  is  the 
preferred  way  of  reducing  the  linearization  error  since  transient  Accuracy  is 
thereby  improved.  Linearization  by  Taylor  series  expansion  in  rime  about  the 
known  time  level  introduces  errors  no  greater  than  those  due  to  the  differencing 
(Refs.  41  and  18),  and  this  approach  has  been  employed  here.  The  formal 
linearization  process  results  in  a  system  of  coupled  equations  in  order  to 
retain  the  requisite  order  of  temporal  accuracy.  The  system  of  coupled 
equations  at  the  implicit  time  level  is  solved  efficiently  using  a  standard 
block  elimination  matrix  inversion  scheme.  In  the  present  problem,  the  strong 
coupling  effects  among  the  governing  equations  dictate  the  use  of  the  block 
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coupled  equation  approach.  However,  weakly  coupled  equations  could  be  solved  in 
a  decoupled  manner  in  order  to  reduce  computer  time  and  storage  requirements.  A 
description  of  the  linearization  technique  and  the  ADI  procedure  is  given  in 
Appendix  B. 

Initial  and  Boundary  Conditions 

The  solution  procedure  requires  that  initial  and  boundary  conditions  be 
specified.  Although  the  time  asymptotic  solution  of  the  Navier-Stokes  equations  is 
independent  of  initial  conditions,  a  prudent  choice  of  initial  conditions  could 
hasten  convergence.  The  initial  conditions  that  were  chosen  for  the  problems  of 
interest  essentially  involved  the  translation  of  the  profiles  of  the  dependent 
variables  at  the  upstream  boundary  throughout  the  computational  domain  with  the 
constraint  that  the  profiles  satisfy  the  boundary  conditions  at  the  surface  of 
the  body. 

The  choice  of  boundary  conditions  is  not  as  arbitrary,  and  depends  on  the 
physics  of  the  problem.  Since  the  external  flow  is  supersonic  at  the  inflow 
boundary,  profiles  of  the  dependent  variables  are  specified  in  accordance  with  an 
inviscid  characteristic  analysis.  At  the  surface  of  the  body  no-slip  and  constant 
wall  enthalpy  boundary  conditions  are  used.  The  latter  condition  was  used  due  to 
its  simplicity,  and  could  have  easily  been  replaced  by  a  derivative  condition 
(heat  transfer  rate)  on  option.  An  additional  condition  is  needed  for  the  remaining 
variable,  the  density.  Both  the  normal  momentum  equation  evaluated  at  the  solid 
boundary  and  the  normal  derivative  of  pressure  have  been  used  successfully.  Along 
the  axis  of  the  boattail,  symmetry  boundary  conditions  are  employed  while  at  the 
downstream  outflow  boundary  second  derivative  extrapolation  conditions  are  applied. 
These  extraneous  conditions  are  used  in  lieu  of  special  differencing  to  eliminate 
the  exterior  points.  In  view  of  the  reasonable  nature  of  the  physical  approximations 
involved  this  is  found  to  be  both  convenient  and  realistic,  and  we  had  no  dif¬ 
ficulty  in  applying  them.  However,  at  the  top  boundary,  where  waves  generated  in 
the  interaction  process  in  the  vicinity  of  the  boattail  must  exit  the  computational 
domain  special  attention  must  be  given  to  the  appropriate  choice  of  boundary 
conditions.  In  order  to  prevent  these  waves  from  reflecting  from  the  top  boundary 
back  into  the  computational  domain  and  perhaps  polluting  the  solution,  a  Mach 
wave  extrapolation  condition  was  developed.  A  description  of  this  boundary 
condition  is  given  in  some  detail  in  the  next  section. 


24 


Mach  Wave  Extrapolation 


12  3 

Consider  a  nonorthogonal  coordinate  system  (y  ,  y  ,  y  )  and  a  Cartesian 

12  3  12  3 

system  (x  ,  x  ,  x  )  with  Cartesian  velocity  components  (u  ,  u  ,  u  ) . 


Through  point  0  we  may  draw  the  streamline  and  Mach  wave  Inclined  at  an 

angle  u  with  respect  to  the  streamline.  The  streamline  is  at  an  angle  0 

3 

measured  with  respect  to  the  Cartesian  coordinate  x  . 


H- o  =  Sln 


(-k) 


(41a) 


e°  * ,an'’  (H 


(41b) 


which  leads  to  a  total  included  angle  with  respect  to  xJ  of 


t  =  do+  H-o 


The  mach  wave  passing  through  point  0  can  be  represented  as  the  vector  s 


T  =  sin  \f/  i.  +  cosv^  i. 
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where  ^  and  ^  are  unit  vectors  in  the  x  and  x  directions  respectively. 

For  any  scalar  function  $  the  gradient  is  given  by 


=  4>„  ?,  + 


(44) 


The  component  of  V$  in  the  s  direction  is  obtained  by  the  dot  product 

T-  V<3>  =  sin^  <£,,  +  cos 


so  that  the  mach  wave  extrapolation  condition,  i.e.  4>  along  a  mach  wave  is 
constant,  becomes 


Sin  +  COS =  0 


(46) 


Relating  the  Cartesian  coordinates  (x*)  to  the  physical  coordinates  (y*) 


we  obtain 


4> 


d<£ 

lx!. 

+  i* 

Jxi 

dy' 

dx‘ 

dy3 

dx' 

(47) 

d$ 

dy' 

dd) 

dy3 

dy* 

dx3 

dy3 

dx3 

Substituting  Eq.  (47)  into  equation  (46)  we  obtain 


d3> 

dy' 


d<P  . 
dy3 


(48) 


where 


a  = 


dy3  dy3 

17  17 

dy1  dy1 

iJf-'+ai* 


_N_ 

D 


(49) 
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Now  a  is  a  nonlinear  function  of  the  velocity  components  u  ,  u  and  h  the  enthalphy 
via  the  evaluation  of  the  angle  ty.  We  must,  therefore,  linearize  this  term  about 
the  known  time  level  n. 


where 


and  similarly 


da  _ 

dll  3 


da 

ah 


n 


Hence 


n+i 

a  - 


(51) 


The  Mach  wave  extrapolation  boundary  condition  is  employed  on  the  first  sweep  of 
the  ADI  procedure  (in  the  x^-  direction),  at  the  *  level.  This  condition  reduces 
to  the  following  equation 


or 


d<P 

ay' 


+  a 


dtp 

a  y3 


0 


a  A <P*  ,  n+.  /_a£\"  .  /  d±_  \"  (52) 

dy1  ^  \  dy3  /  \  dy 1  / 
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using  the  linearized  form  of  a 


we  obtain 


n+1 


(53) 


1  3 

where  <j>  is  either  u  ,  u  or  h. 

One  could  lag  the  Mach  angle  (evaluate  it  at  the  n  time  level)  by  setting 
fln  ■  0.  For  our  case  this  simplification  was  used  since  we  were  interested 
only  in  a  steady  solution,  and  the  flow  near  the  outer  boundary  was  well 
behaved . 


DISCUSSION  OF  RESULTS 


Several  different  aft  end  flow  fields  are  discussed.  In  order  to  verify 
the  code,  i.e.  L-shaped  domain  and  comer  point  treatment  as  well  as  the 
adaptive  grid  capability  the  laminar  supersonic  two-dimensional  right  angle 
backstep  was  considered.  Turbulence  modeling  effects  were  investigated  in  the 
supersonic  axisymmetric  test  case  of  Badranarayan  (Ref.  57).  In  both  these  cases 
the  expected  qualitative  behavior  was  obtained.  However,  the  results  will  not  be 
discussed  in  any  great  detail  except  where  relevant  to  the  studies  they  were 
specifically  designed  to  address.  The  cases  that  are  considered  in  detail  are 
the  flow  over  the  AGARD  10°  nozzle  boattail  at  M  =  1.5  with  and  without  let 
exhaust.  These  cases  correspond  to  the  experiments  of  Galigher,  et  al  (Ref.  5). 

In  addition,  the  two-dimensional  supersonic  flow  over  a  rearward  facing  step  and 
the  subsequent  reattachment  of  the  free  shear  layer  on  a  20°  inclined  ramp  will 
also  be  described.  This  case,  which  corresponds  to  the  data  of  Settles,  et  al 
(Ref.  58)  was  submitted  by  SRA  to  the  1981  AFOSR  -HTTM-  Stanford  Conference  on 
Complex  Turbulent  Flows. 

Mesh  Movement  -  Adaptive  Grid 

In  the  solution  of  problems  where  regions  of  large  gradients  occur,  lack  of 
mesh  resolution  can  lead  to  oscillations,  cell  Reynolds  problems  or  give  misleading 
results.  Whenever  these  regions  can  be  isolated  and  identified  then  standard  mesh 
refinement  techniques  such  as  those  discussed  in  this  report  can  be  implemented. 
However,  in  many  cases  these  regions  either  cannot  be  identified  a  priori  or  they 
change  in  time  as  the  solution  procedure  progresses.  There  is,  therefore,  great 
interest  in  being  able  to  adapt  the  grid  to  the  particular  problem.  Two  aspects 
of  the  grid  adaptor  procedure  must  be  considered.  The  first  is  the  actual  process 
of  moving  the  grid  and  the  second,  the  far  more  difficult,  is  the  choice  of 
criteria  by  which  the  grid  moves.  The  latter  point  becomes  even  more  difficult 
if  we  wish  to  automatically  adapt  the  grid  by  monitoring  particular  functions  or 
derivatives,  since  we  do  not  wish  to  mislead  the  geometry  into  clustering  grid 
points  in  inappropriate  regions  or  by  moving  it  too  rapidly. 

The  first  problem  of  actual  grid  movement  is  accomplished  by  using  the  time- 
dependent  geometry  form  of  Navier-Stokes  equations  that  has  been  applied  by  Gibeling 
and  McDonald  Ref.  (  59  )  to  interior  ballistics.  The  metric  data  is  cast  into  a 
form  that  it  is  a  function  of  time.  Hence,  the  governing  equations  contain  terms 
that  account  for  geometry  movement  (they  appear  as  convection  terms) .  This 
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procedure  eliminates  any  need  to  interpolate  data  from  one  grid  onto  another 
across  a  time  step.  An  alternate  approach  would  entail  interpolating  onto  a  new 
grid  at  the  end  of  a  time  step,  and  then  marching  the  solution  from  there.  The 
main  advantage  of  this  procedure  is  that  it  simplifies  the  evaluation  of  geometric 
germs  since  it  is  now  independent  of  time.  However,  it  has  the  disadvantage, 
as  previously  mentioned,  of  requiring  the  interpolation  of  the  entire  flow  field 
onto  the  new  grid  whenever  the  grid  is  moved.  Since  the  time  dependent  geometry 
option  allows  for  geometries  that  can  change  in  time,  i.e.  expanding  domains,  the 
additional  generality  influenced  the  choice  of  method.  However,  the  alternate 
approach  could  be  useful  for  certain  applications  and  deserves  further  consideration. 

The  second  problem  of  choosing  the  criteria  by  which  the  mesh  is  modified 
in  time  is  now  considered.  Our  aim  is  to  redistribute  the  number  of  grid  points 
so  that  they  are  concentrated  in  the  region  of  the  free  shear  layer.  In  order 
to  locate  the  shear  layer,  the  region  of  large  gradients,  an  optimization  routine 
developed  by  Levy  (Ref.  60)  was  adapted  for  the  MINT  code.  We  assume  that  the 
shear  layer  begins  at  the  corner  point  (cf.  Fig.  3)  and  can  be  represented  by  an 
analytic  function.  In  this  case,  it  was  chosen  to  be  a  cosine  function,  i.e. 

x*  xQ  +  A  ( I  -  cosXCz  -  ZQ))  X(z— Zg)£ll  (54) 

X»Xq+2A  X  (z-  Zq)  >  n 

where  x  and  z  are  the  coordinates  of  the  corner  and  the  two  underdetermined 
o  o 

coefficients  A  and  X,  the  amplitude  and  wave  length  respectively  are  functions  of 
time.  The  cosine  function  was  chosen  due  to  its  simplicity,  but  any  other 
convenient  function  could  have  been  used .  In  the  vicinity  of  the  shear  layer  the 
flow  field  is  sampled  for  the  magnitude  of  the  first  derivative  of  the  dependent 
variables.  The  data  is  then  filtered  and  an  "optimum"  curve  is  derived  (e.g.  A  and 
X  are  chosen)  that  best  fits  the  data  and  has  the  functional  form  given  by 
equation  (54). 

Assume  for  the  moment  that  the  original  mesh  was  orthogonal  and  that  a  Robert's 
type  transformation  was  employed  that  was  centered  about  the  line  x  *  xq,  where  xq 
is  the  location  of  the  corner  (cf.  Fig.  3).  The  shear  layer  will  in  general  depart 
from  the  original  horizontal  line  x  *  xq.  Since  we  wish  to  cluster  points  in  the 
vicinity  of  the  shear  line,  we  center  the  grid  transformation  about  this  new  computed 
curve.  Note  that  the  resulting  grid  will  no  longer  be  orthogonal  in  the  region  of 
the  shear  layer.  Both  levels  of  geometry,  the  original  orthogonal  grid  at  the  n 
time  level  and  the  new  nonorthogonal  grid  at  the  (n  +  l)st  time  level  are  employed  to 
advance  the  solution  one  step  in  time.  At  this  point,  the  data  can  be  sampled  again 
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to  determine  a  new  location  of  the  shear  layer.  The  process  is  repeated  at  some  pre¬ 
determined  frequency  (in  time)  or  until  the  shear  layer  position  does  not  change 
within  some  preset  tolerance. 

This  adaptive  grid  option  was  examined  for  the  laminar  two-dimensional  flow 
over  a  back  step  at  a  Mach  number  of  =  3  and  Reynolds  number  based  on  step 

height  of  Re^  =  400.  A  velocity  vector  plot  of  the  flow  field  is  shown  in  Fig.  4. 

Two  tests  were  conducted.  In  the  first,  the  grid  was  forced  to  move  through 
half  a  cycle  which  corresponds  to  the  maximum  deviation  from  the  initial  state  of 
a  horizontal  line.  The  purpose  of  this  test  was  to  verify  that  the  solution  is 
well  behaved  during  the  mesh  movement,  and  to  determine  how  sensitive  the  solution 
was  to  the  mesh  clustering  about  the  shear  layer.  Figure  5  shows  the  velocity 
profiles  at  two  locations  downstream  of  the  base,  in  the  recirculation  zone  and 
in  the  wake  before  and  after  the  mesh  movement.  The  results  compare  very  well 
with  each  other  indicating  that  the  method  does  not  introduce  errors  due  to  grid 
movement  and  that  in  the  case  of  interest  the  solution  is  not  very  sensitive  to 
the  location  of  the  shear  layer.  It  is  expected,  however,  that  for  other  cases  in 
particular  for  turbulent  flows  greater  mesh  sensitivity  would  be  encountered. 

The  second  test  case  was  designed  to  allow  the  grid  to  move  under  the  control  of  the 
optimization  curve  fit  procedure  described  previously.  In  Fig.  6  the  velocity  pro¬ 
files  in  the  wake  and  recirculation  zone  are  compared  again  prior  to  and  after  mesh 
movement.  The  results  compare  well  with  one  another.  In  the  case  considered,  the 
original  orthogonal  grid  had  adequate  resolution  in  the  shear  layer  and,  therefore, 
the  use  of  an  adaptive  grid  did  not  have  a  substantial  effect  on  the  results.  In 
other  problems  being  considered  at  SRA  this  is  not  the  case.  For  instance  in  the 
normal  shock  calculations  (ARO  Contract  DAAG29-80-C-082)  the  adaptive  grid  procedure 
aided  in  stabilizing  the  moving  shock  and  resolving  the  details  of  its  structure. 

Although  for  the  cases  considered  in  this  report  (e.g.  boattail  geometry)  the 
adaptive  grid  option  would  have  been  useful,  we  did  not  proceed  with  its  implementa¬ 
tion  at  the  present  time.  A  primary  reason  was  that  the  hyperbolic  function  grid 
transformations  that  were  employed  initially  in  order  to  cluster  points  within  the 
shear  layer  did  not  have  the  desired  grid  control  under  mesh  movement.  This  is 
in  part  a  consequence  of  the  physics  of  the  problem  which  leads  to  a  "shear  layer 
curve"  that  is  a  function  of  x  and  z  as  compared  to  the  normal  shock  which  is 
primarily  a  function  of  z  alone.  For  the  normal  shock  case  adequate  grid  control 
was  obtained.  Recently,  new  grid  transformations  have  been  developed  that  alleviate 
this  problem  and  allows  for  better  grid  control.  They  have  not  as  yet  been  used  in 
adaptive  grid  calculations. 
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Nozzle  Boattail 


Experimental  data  is  very  limited  for  the  cases  of  interest;  the 
supersonic  axisymmetric  turbulent  flow  over  a  boattail  geometry  with  and 
without  jet  exhaust.  Either  the  geometry  is  inappropriate,  i.e.  2-D  versus 
axisymmetric,  or  the  flow  is  subsonic  or  insufficient  data  is  given.  The  case 
chosen  with  which  to  compare  our  calculation  was  that  of  Galigher,  et  al  (Ref.  5). 
This  case  has  two  main  shortcomings,  the  data  is  incomplete,  there  are  only 
pressure  measurements  for  the  shoulder  and  no  measurments  at  all  in  the  interaction 
zone  downstream  of  the  base,  and  there  are  three  dimensionality  effects  due  to  the 
model  supports.  Although,  at  the  outset  of  our  calculation  these  shortcomings 
were  evident,  it  was  felt  that  the  data  for  the  pressure  signature  on  the  shoulder 
should  be  reliable,  and  of  sufficient  interest  to  warrant  their  consideration. 
Furthermore,  the  jet  on  case  was  computed  by  Mikhail,  et  al  (Ref.  22),  and  this 
seemed  to  afford  us  a  means  of  comparing  our  calculations  with  other  results. 
Unfortunately,  the  lack  of  data  was  a  far  greater  detriment  than  expected. 

Since  our  results  indicate  that  the  pressure  signature  on  the  shoulder  is 
sensitive  to  the  turbulence  model,  additional  data  on  the  shoulder,  i.e.  velocity 
profiles  as  well  as  in  the  jet  and  wake  could  have  been  useful  in  selecting  the 
appropriate  length  scales.  Although  the  calculations  showed  less  sensitivity  to 
upstream  turbulent  boundary  layer  parameters,  additional  data  by  the  experimentalists 
could  have  been  of  value.  Originally,  we  had  intended  to  compare  our  calculations 
with  those  of  Mikhail,  et  al  (Ref.  22).  However,  due  to  significant  differences  in 
the  geometry  adopted  and  substantial  differences  in  turbulence  modeling  as  compared 
to  that  adopted  here,  it  was  decided  that  a  comparison  with  the  calculation  of  Mikhal 
could  be  misleading  and  this  is  deferred  until  some  of  the  questions  that  are  raised 
later  in  this  report  can  be  resolved. 

In  this  section  we  describe  the  construction  of  the  coordinate  system  for 
the  boattail  geometry,  the  boundary  conditions  and  turbulence  model  employed  and 
the  results  of  our  calculations  for  the  jet  on  and  jet  off  cases. 

Figure  7  shows  the  AGARD  10®  nozzle  boattail  configuration  with  the 
coordinate  information  given  in  tabular  form.  In  order  to  simplify  subsequent 
computations,  the  geometric  data  was  represented  by  a  least  square  cubic  curve 
fit  given  as 

r  ■  4.93  -  .022948X  -  .00302X2  -  .0003747X3  (55) 
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The  maximum  error  between  the  analytic  representation  and  the  actual  boattail 
shape  was  less  than  1.5%. 

The  computational  domain  is  shown  in  Fig.  8,  and  includes  moving  counter¬ 
clockwise  the  upstream  boundary,  boattail  shoulder,  boattail  base  and  jet  exit, 
symmetry  line,  downstream  boundary  and  the  outer  boundary.  Note  that  the  region 
external  to  the  boattail  body  is  L-shaped.  The  upstream  boundary  was  chosen  at 
5  baseheights  upstream  of  the  base,  the  downstream  boundary  at  8  base  heights 
downstream  of  the  base  and  the  outer  boundary  at  4  base  heights  above  the  symmetry 
line. 

The  first  task  was  to  construct  a  coordinate  system  with  a  sufficient 
distribution  of  grid  points  to  adequately  resolve  the  flow  field  within  the 
computational  domain,  while  satisfying  the  constraints  that  all  boundaries  must 
lie  on  coordinate  lines  and  the  base  shoulder  interaction  remain  sharp.  As 
discussed  in  an  earlier  section,  this  leads  to  a  nonorthogonal  grid  distribution. 
In  the  following,  the  grid  construction  is  outlined. 

Consider  the  rectangle  enclosing  the  boattail.  It  includes  as  sides 
the  upstream  and  downstream  boundaries  and  the  outer  boundary  and  symmetry 
lines.  The  z  =  constant  lines  (vertical  lines)  are  constructed  first. 


Since  these  lines  do  not  vary  with  x  a  single  distribution  is  required  which  j 

assures  that  the  upstream  and  downstream  boundaries  lie  on  grid  points  and  that  | 

the  base,  which  is  in  the  interior  also  lies  on  a  grid  point.  We  employ  an 
analytical  transformation  which  satisfies  the  above  conditions  and  allows  for 
grid  spacing  control  to  the  requisite  degree  of  resolution  in  the  vicinity  of 
the  base.  Refer  to  Fig.  9  for  the  z  distribution.  In  the  case  we  considered,  j 

there  were  41  grid  points  in  the  z  direction  with  the  base  falling  on  the  20th  j 

grid  point.  j 

Having  obtained  the  z  distribution  we  are  in  a  position  to  construct  the  j 

x  distribution.  The  constraints  for  this  portion  of  the  grid  generation  process 
are  that  the  outer  and  symmetry  boundaries  lie  on  grid  points.  In  addition,  j 

the  boat-tail  surface  must  also  be  a  coordinate  line.  Since  the  boattail  j 

surface  terminates  at  the  base,  we  must  analytically  extend  it  in  order  that  the  metric 
data  be  continuous  at  the  corner  point.  This  was  accomplished  by  passing  a  second 
order  polynomial  through  the  corner  point,  which  matched  the  slope  of  the  boattail 
at  the  corner  and  terminated  at  some  location  downstream  of  the  base  with  a  zero 
slope.  For  each  z  location,  the  upper  boundary,  the  boattail  surface  and  its 
extension  and  the  symmetry  line  were  required  to  lie  on  grid  lines  with  the 
additional  constraint  that  the  grid  spacing  near  the  boattail  surface  be  chosen 
to  adequately  resolve  the  turbulent  sublayer.  The  analytic  transformation  used  in 
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obtaining  the  z  distribution  was  also  used  here  and  obtained  the  grid  lines  shown 
in  Fig.  9.  Note  that  in  contrast  to  the  z  distribution  which  was  not  a  function  of  x 
and  was  obtained  in  one  pass,  the  x  distribution  was  repeated  at  every  z  station. 

This  is  a  direct  consequence  of  the  nonorthogonality  of  the  grid,  i.e.  the  x 
coordinate  lines  are  functions  of  z.  In  the  case  under  consideration  41  grid 
points  were  used  in  the  x  direction  with  the  boattail  surface  lying  on  the 
20th  grid  point.  Since  the  points  within  the  nozzle  boattail  do  not  lie  within 
the  computational  domain,  a  total  of  1320  grid  points  including  boundaries  were 
used  for  the  computation.  This  is  a  rather  meager  number  of  grid  points,  and  large 
stretches  had  to  be  used  in  order  to  satisfy  the  constraints  that  were  placed  on  the 
transformation.  More  grid  points  could  have  eased  the  severity  of  the  stretches 
while  allowing  better  resolution.  Whether  there  was  adequate  resolution  in  the 
interaction  zone  is  still  open  to  question.  The  jet  on  results  indicate  that 
immediately  upstream  of  the  base  there  was  insufficient  resolution.  This  question 
will  be  addressed  in  greater  detail  below. 

In  order  to  solve  the  Navier-Stokes  equtions,  appropriate  boundary 
conditions  must  be  specified.  These  are  given  below. 

1.  Upstream  inflow  boundary 

Profiles  specified  for  u,  w,  h,  p 

2.  Boattail  surface  and  base 

u  =  w  =  0 

T  =  Tw 


3.  Symmetry  line 

u  *  0 

*z.  .  m  iP  .  o 

3n  3n  3n 

4.  Downstream  outflow  boundary 

2  2  2  2 

3  Si  m  3  w  .  3  n  _  3  P  _  Q 

2  2  2  2 
3z  3z  3z  3z 

5.  Top  outflow  boundary 

Mach  wave  extrapolation  for  u,  w,  h,  p 

6.  Jet  at  base 

u  -  0 

w,  h,  p  or  P  profiles  specified 
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The  wall  temperature  was  taken  as  Tw  =  580°R.  In  order  to  obtain  upstream 
turbulent  profiles,  the  boundary  layer  thickness  is  needed.  A  value  of  6  =  1.12 
inches  was  chosen  which  corresponds  to  the  estimate  of  Mikhail,  et  al  (Ref.  22). 
Insufficient  information  is  given  in  the  experiment  to  obtain  a  precise  value. 

A  better  approximation  for  the  upstream  profiles  and  6  could  be  obtained  by 
exercising  a  marching  procedure,  e.g.  the  PEPSI-S  Code  (Ref.  61)  and  marching  from 
the  nose  of  the  one  cylinder  until  the  upstream  boundary  of  the  Navier-Stokes 
computational  domain.  Numerical  experiments  indicated  that  decreasing  6  did  not 
have  an  appreciable  effect  on  the  pressure  signature  on  the  shoulder. 

Employing  the  procedure  of  Maise  and  McDonald  (Ref.  62)  the  appropriate 
constants  for  the  Coles  law  profile  can  be  obtained,  and  streamwise  velocity 
profile  computed.  Assuming  Hq  =  constant,  u  *  0,  and  that  the  pressure  is 
constant  across  the  boundary  layer  and  equal  to  the  free  stream  value,  the  density 
and  enthalpy  profiles  can  be  obtained.  Since  these  profiles  are  not  precisely 
consistent  with  the  interior  solutions,  a  procedure  was  implemented  in  order  to 
eliminate  their  effect  on  the  Navier-Stokes  solution.  A  preliminary  calculation 
was  performed,  beginning  upstream  of  the  Navier-Stokes  region,  at  z/h  <*  -8  and  term¬ 
inating  at  z/h  =-3.  This  computational  domain  starts  in  uniform  flow  and  has  no 
knowledge  of  the  base  and  hence  will  not  be  influenced  by  upstream  effects.  The 
Navier-Stokes  equations  are  solved  within  this  domain,  and  the  profiles  that  are 
obtained  at  z/h  =  -5  then  become  the  upstream  boundary  conditions  for  the  base 
region  domain.  The  pressure  coefficient  obtained  for  this  calculation  is  presented 


in  Fig.  10.  Also,  on  this  figure  are  the  c  curves  for  the  jet  on  condition.  Note 
the  three  dimensionality  of  the  experiment  in  the  variation  of  the  c^  curves  for 
the  upper  and  lower  surfaces.  The  computed  c  curve  which  lies  between  the  two 
experimental  data  curves  is  in  fairly  good  agreement  with  the  data. 


A  turbulence  model  must  now  be  specified  in  order  to  obtain  closure.  An 
algebraic  mixing  length  model,  the  one-equation  k-£.  model  and  the  two-equation  k-e 
model  were  experimented  with  for  a  number  of  test  cases.  In  our  calculations  of  the 
supersonic  axisymmetric  backstep  problem  of  Badranarynan  (Ref.  57),  a  comparison  of 
the  mixing  length  and  k-X.  models  was  made.  Both  models  are  dependent  on  evaluating 
an  appropriate  length  scale.  On  the  surface  of  the  body,  this  is  relatively  easy. 
However,  in  the  wake  it  Is  much  more  difficult  since  we  must  detect  the  "edge"  of 
the  shear  layer  which  is  not  well  defined.  However,  the  results  obtained  by  both 
methods  in  the  wake  of  the  body  were  very  close  to  one  another. 
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During  the  calculation  of  the  Stanford  test  case  (see  following  section) 
the  two-equation  k-e  model  was  exercized.  Difficulties  were  encountered  immediately 
downstream  of  the  corner  in  the  recirculation  zone.  At  least  some  of  these  dif¬ 
ficulties  were  traced  to  the  low  Reynolds  number  correction  terms  that  were  devised 
to  be  used  near  solid  boundaries,  which  apparently  give  rise  to  non-physical  be¬ 
havior  in  the  free  shear  layer.  In  view  of  the  lack  of  confidence  in  the  k-e  model 
with  no-slip  boundary  conditions  and  the  absence  of  any  substantial  improvements  in 
using  the  computationally  more  expensive  k-Jl  model,  the  algebraic  mixing  length 
model  was  used  to  obtain  the  reported  results. 

Figures  11  -  16  show  computer  plots  of  the  results  of  the  calculation. 

In  Fig.  11  the  Mach  number  contours  for  the  entire  computational  domain  are 
presented.  The  interaction  of  the  expansion  wave  on  the  shoulder  with  the 
separation  compression  wave  is  evident.  Furthermore,  the  closing  of  the  Mach 
number  contours  incidate  the  presence  and  extent  of  the  recirculation  zone. 

A  blown  up  view  of  the  coordinate  system  in  the  vicinity  of  the  base  is  presented 
in  Fig.  12.  Several  streamwise  coordinate  lines  near  the  shoulder  of  the 
shoulder  of  the  boattail  have  been  eliminated  for  clarity.  The  pressure  contours 
in  this  region  are  shown  in  Fig.  13.  Note  that  the  pressure  minimum  is  in  the 
center  of  the  recirculation  vortex  and  that  the  flow  recompresses  downstream 
of  this  location  in  the  near  wake.  Further  downstream  (not  shown  in  this  figure) 
in  the  far  wake,  the  flow  reexpands.  The  next  two  figures  14  and  15  show  the 
velocity  vector  plot  for  the  full  domain  and  the  near  wake  region  respectively. 

The  turning  of  the  shear  layer  and  the  resulting  closed  recirculation  zone  is 
clearly  seen.  These  results  are  in  keeping  with  the  expected  qualitative 
behavior. 

As  mentioned  earlier,  the  sole  data  presented  in  Ref.  5  is  for  the  pressure 
coefficient,  cp,  on  the  shoulder  of  the  boattail.  In  Fig.  16  we  present  a 
comparison  of  our  calculations  with  the  data  of  Ref.  5.  From  1.5  to  8  base 
heights  upstream  of  the  base,  the  data  for  the  upper  and  lower  surfaces  are 
shown  with  the  spread  in  cp  values  being  indicative  of  the  three  dimensionality 
of  the  flow.  For  the  section  of  the  boattail  in  the  vicinity  of  the  comer  only 
the  top  surface  data  is  shown  since  we  did  not  have  available  the  bottom  surface 
data.  In  the  region  from  3  to  5  base  heights  upstream  the  computed  cp  curve 
is  indistinguishable  from  the  curve  presented  in  Fig.  10  which  was  used  to  obtain 
the  upstream  boundary  conditions  and  confirms  that  the  location  of  the  upstream 
boundary  placement  is  not  influenced  by  the  interaction  process. 


As  compared  to  the  experimental  data,  the  computations  show  a  pressure  rise 
that  is  delayed  and  thus  leads  to  a  smaller  region  of  upstream  influence.  Although 
the  pressure  rise  has  the  same  slope  as  the  data,  there  is  no  plateau  or  leveling 
out  of  the  pressure  as  the  corner  is  approached. 

Several  factors  could  have  contributed  to  the  quantitative  lack  of  agreement. 

is  the  turbulence  modeling  in  the  wake  and  on  the  shoulder.  The  turbulence  model 

employed  was  an  algebraic  mixing  length  model.  The  mixing  length  on  the  shoulder 

was  computed  by  Eq.  (27).  In  the  free  shear  layer  and  wake  the  mixing  length  was 

nominally  computed  from  a  linear  expression  relating  at  the  corner  and  J l 

at  some  location  downstream  in  the  wake.  The  value  of  SL  at  the  corner  was  the 

00 

value  obtained  at  the  point  of  minimum  pressure  on  the  boattail  shoulder  and 
which  was  frozen  at  that  value  over  the  remainder  of  the  boattail  surface. 

In  order  to  assess  the  effect  of  the  length  scale  in  the  wake  region  on  the 
boattail  corner  pressure,  the  value  of  l  was  varied.  This  is  shown  in 
Fig.  17.  As  Z*  is  diminished,  the  effective  JL^  in  the  recirculation  is 
increased  resulting  in  a  greater  turbulent  viscosity  and,  hence,  a  lower 
corner  pressure.  The  c  results  presented  in  Fig.  16  are  for  case  4. 

The  effect  of  the  i  distribution  on  the  boattail  surface  was  not 
00 

evaluated  in  such  detail,  but  it  would  appear  from  our  computations  that  it 

has  a  less  significant  effect  than  the  distribution  in  the  wake.  However, 

the  inherent  assumptions  of  the  mixing  length,  i.e.  the  turbulent  viscosity  is 

proportioned  to  the  velocity  gradients  may  be  suspect  in  the  separated 

region  on  the  boattail.  Further  detailed  investigation  of  the  one  and  two-equation 

turbulence  models  would  be  of  value  in  resolving  this  issue. 

It  must  also  be  pointed  out  that  the  jet  off  computation  did  not  correspond 
precisely  with  the  actual  experiment.  Apparently,  the  jet  off  experiment  was 
run  with  the  nozzle  open.  This  would  allow  fluid  to  enter  the  nozzle  and  thus 
give  rise  to  a  recirculation  zone  that  was  different  than  one  that  would  exist  if 
nozzle  plane  was  a  solid  surface.  The  flow  structure  in  the  vicinity  of  the  base 
would  necessarily  have  an  effect  on  the  pressure  at  the  corner. 

To  calculate  the  flow  field  with  the  nozzle  open  would  require  that  the 
computational  domain  include  a  region  inside  the  nozzle  upstream  of  the  exit  plane. 
Although  the  computer  code  is  designed  to  handle  such  geometries  at  this  stage  of 
the  model  development  it  was  felt  that  the  added  complexity  was  not  warranted. 

Thus,  in  the  present  computations  a  solid  wall  on  the  nozzle  exit  plane  was  assumed. 
In  comparing  the  computations  with  data,  it  would  appear  that  assuming  a  solid 
wall  could  have  a  significant  effect  on  the  boattail  pressure  signature.  It  is 
recommended  that  in  any  study  aimed  at  resolving  the  discrepencies  in  the  c^ 
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distribution  additional  effort  should  include  expanding  the  computational  domain 
into  the  interior  of  the  nozzle. 


Jet  On  Case 

The  next  case  we  considered  was  a  nozzle  boattail  with  jet  exhaust  at  a 

nozzle  pressure  ratio  (NPR)  of  P.  „/P  =7.09.  P  is  the  jet  static  pressure 

jet  oM  jet 

and  Pq  is  the  free  stream  stagnation  pressure.  The  other  pertinent  data  are: 


T 


jet 


450°R 


M 


'jet 


1.08 


Nozzle  exit  radius  R  =  1.991  inches 


Since  the  jet  is  supersonic, function  data  in  the  form  of  profiles  of  the 
dependent  variables  can  be  applied  as  boundary  conditions.  Note  that  the  jet 
does  not  span  the  entire  base  area  and  this  results  in  split  boundary  conditions 
on  the  base  that  changes  abruptly  at  the  jet/wall  juncture.  No  computational 
problems  were  observed  as  a  result  of  applying  such  boundary  conditions. 

The  jet  issuing  from  the  nozzle  is  assumed  to  have  a  streamwise  velocity 
distribution  given  by 


W{  r)  =  W.  0<  r  <  R* 

/  □  _r  \  1/7  * 

W(r)  =  Wj  R  <  r  <R 


(56) 


where  W.  is  the  jet  centerline  velocity,  6  =  ,1R  and  R  =  .9R,  while  the  normal 
velocity  is  set  to  zero.  In  addition  the  stagnation  temperature  of  the  jet  is 
assumed  constant  as  is  the  static  pressure.  This  model  was  used  to  obtain  the 
appropriate  density  and  static  enthalpy  profiles.  The  same  grid  distribution  was 
used  for  both  the  jet  off  and  jet  on  cases.  At  the  base,  the  jet  exit  area 
encompassed  grid  points  1-11  while  the  solid  wall  was  located  on  grid  points  12  -  20. 

The  results  of  the  calculations  are  shown  in  Figs.  18  -  22.  In  Fig.  18 
the  Mach  number  contours  are  presented  for  entire  computational  domain.  The 
pattern  is  very  different  from  the  jet  off  case.  With  the  jet  on  there  is  very 
little  turning  and  in  the  near  wake  there  is  a  rapid  expansion.  This  is  more 
clearly  evident  in  Fig.  19  where  the  pressure  contours  are  shown  in  the  vicinity 
of  the  base.  In  Figs.  20  and  21  the  velocity  vector  plots  for  the  entire 
computational  domain  and  for  the  near  wake  region  are  presented.  The  recirculation 
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zone  on  the  shoulder,  upstream  of  the  corner,  can  be  seen  indicating  the  relatively 
large  region  of  upstream  influence.  These  results  are  in  keeping  with  the  expected 
physical  behavior  of  the  jet  interaction. 

The  computed  cp  distribution  on  the  shoulder  is  shown  in  Fig.  22  where  it  is 
compared  to  the  data  in  Ref.  5.  The  pressure  rise  in  the  vicinity  of  the  base  is 
somewhat  delayed  when  compared  to  the  data,  and  overshoots  the  experimental  peak  cp. 
Also,  near  the  corner  oscillations  in  cp  are  noted.  The  oscillations  are  believed 
grid  dependent  and  could  be  eliminated  by  resolving  the  region  near  the  corner  with 
a  finer  grid.  The  reason  for  the  overshoot  and  delay  in  the  pressure  rise  is  still 
unresolved,  but  as  in  the  jet  off  case,  it  is  believed  that  turbulence  modeling  plays 
a  significant  role. 

Rearward  Facing  Step 

The  solution  procedure  described  in  the  previous  section  is  also  used  to 
compute  the  two-dimensional  supersonic  turbulent  flow  over  a  rearward  facing  step 
with  shear  layer  reattachment  on  an  inclined  ramp.  The  numerical  results  that  are 
compared  with  the  data  of  Settles,  et  al  (Ref.  57)  were  submitted  by  SRA,  Inc. 
to  the  "1981  AFOSR  -HTTM-  Stanford  Conference  on  Complex  Turbulent  Flows". 

The  interesting  feature  of  this  problem  is  that  the  process  of  reattachment  can  be 
studied  without  complicating  the  upstream,  incoming  flow,  with  shock  and  expansion 
wave  effects.  These  effects  are  present  when  we  consider  a  simple  back  step  flow. 

In  that  case  the  expansion  wave  eminating  from  the  vicinity  of  the  corner  and  the 
recompression  shock  influence  the  shear  layer  development  and  complicate  the 
initial  conditions  of  reattachment.  By  carefully  adjusting  the  experiment 
Settles,  et  al  (Ref.  58)  were  able  to  isolate  the  reattachment  process  by  allowing 
the  free  shear  layer  to  leave  the  back  step  parallel  to  the  wall  with  no  wave 
disturbance. 

The  experimental  model  and  a  schematic  of  the  flow  field  is  shown  in  Fig.  23 

with  the  appropriate  dimensions.  The  flow  field  contains  two  corners,  the  back  step 

and  the  cavity  ramp  juncture.  The  former  was  kept  sharp  since  the  geometry  at  that 

location  has  a  strong  influence  on  the  flow  structure  while  the  ramp  cavity  juncture 

was  rounded  with  a  small  circular  fillet.  The  rounding  of  the  cavity  ramp  juncture 

is  not  expected  to  be  an  appreciable  effect  on  the  structure  of  the  redeveloping 

flow.  The  construction  of  the  geometry  deserves  special  attention.  Since  the 

incoming  flow  is  parallel  to  the  upper  surface  of  the  back  step,  it  is  desirable 
3 

to  construct  the  x  coordinate  lines  parallel  to  this  surface.  Similarly,  since 
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the  outgoing  flow  will  eventually  align  itself  with  the  ramp,  it  is  also  desirable 
3 

to  construct  the  x  coordinate  lines  parallel  to  the  ramp  near  the  exit  of  the 

3 

computational  domain.  Within  the  computational  domain  the  x  coordinate  lines 

must  then  vary  smoothly  between  these  two  extremes,  the  x^  coordinate  lines 

3 

by  the  same  token  should  be  orthogonal  to  the  x  coordinate  lines  at  the  entrance 
and  exit  of  the  computational  domain.  However,  within  the  computational  domain 
the  coordinates  become  nonorthogonal  (cf.  Figs.  24  and  25).  The  actual  construction 
of  the  coordinate  system  is  briefly  described  below. 

Consider  the  computational  domain  (cf.  Fig.  24)  which  is  set  off  by  the  bold 
lines.  Moving  counterclockwise  we  start  at  the  upstream  boundary,  and  proceed 
over  the  back  step  surface,  the  base  and  the  cavity  floor.  This  is  followed  by 
a  20°  inclined  ramp,  the  exit  plane  that  is  perpendicular  to  the  ramp  and  finally 
the  outer  boundary  which  begins  parallel  to  the  ramp  and  ends  up  parallel  to  the 
back  step  surface.  Note  that  the  cavity  ramp  juncture  as  well  as  the  outer 
boundary  are  rounded  with  circular  fillets  so  that  the  metric  data  in  these 
locations  would  remain  smooth.  In  order  to  construct  the  desired  coordinate 
system,  we  focus  on  the  trapezoid  ABCD  that  encompasses  the  computational  domain. 
First  consider  the  coordinate  lines  which  intersect  the  upper  and  lower 
boundaries.  The  constraints  which  are  imposed  are  that  the  inlet  AED,  the  exit 
BGC  and  the  line  containing  the  base  all  lie  on  coordinate  lines.  Furthermore,  with 
the  conditions  that  the  grid  spacing  is  prescribed  these  lines  and  that  the 
line  adjacent  to  the  exit  is  orthogonal  to  the  ramp  we  can  construct  the  required 
grid  distribution  by  employing  a  hyperbolic  function  transformation.  This  is 
accomplished  by  distributing  points  both  along  the  top  and  bottom  sides  of  the 
trapezoid  (lines  AB  and  CD)  and  connecting  the  corresponding  points  by  straight 
lines  (cf.  figures  24  and  25). 

In  order  to  obtain  the  other  family  of  coordinates  we  use  the  same  approach 
as  for  the  boattail  and  construct  the  curve  EFG  which  has  zero  slope  at  the  corner 
F  and  is  parallel  to  the  ramp  at  point  G.  Employing  the  procedure  used  in  the 
construction  of  the  streamwise  cooridnate  lines  for  the  boattail  geometry,  the 
coordinate  system  shown  in  Fig.  25  is  obtained. 

The  flow  conditions  for  the  problem  are 

M  -  2.92 
00 

Re  *  6.7  x  10^/meter 

Tw  -  265°k 
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The  governing  equations,  continuity,  two  Cartesian  momentum,  energy  and 
state  are  solved  coupled,  with  the  following  boundary  conditions: 

1.  Upstream  inflow  boundary 

Profiles  specified  for  u,  w,  h,  p 

2.  Backstep  surface,  base,  cavity  and  ramp 

u  =  w  =  0 
T  =  Tw 


3.  Downstream  outflow  boundary 

2  2  2  2 

9  u  =  JTw  _  JTh  _  2_P. 

2  2  2  ~  2 
9z  9z  9z  9z 

4.  Top  outflow  boundary 

Mach  wave  extrapolation  for  u,  w,  h,  p 

The  turbulence  model  used  here  was  also  an  algebraic  mixing  length  model  with 
linearly  increased  in  the  wake  until  it  reached  a  value  of  =  .09  <$a  the 
average  6  on  the  ramp  which  corresponds  to  Z  *  .039. 

It  has  been  found  necessary  to  add  an  artificial  dissipation  term  in  the 
numerical  algorithm  in  order  to  suppress  wiggles  that  could  appear  in  regions  of 
large  cell  Reynolds  number.  However,  adding  too  much  dissipation  could  smear 
out  shock  waves.  In  Ref.  (63)  it  is  shown  that  by  appropriately  dialing  down  the 
dissipation  we  can  sharpen  the  shock  waves.  Such  a  method  has  been  employed  for 
the  results  presented  here. 

In  Fig.  (26)  the  pressure  contours  are  shown.  The  pressure  disturbance 
produced  at  the  corner  is  inclined  at  the  local  Mach  angle,  while  the  compression 
wave  that  is  generated  in  the  reattachment  process  coalesces  into  a  shock.  The 
velocity  vector  plot  presented  in  Fig.  (27)  confirms  that  the  shear  layer  is 
parallel  to  the  cavity  wall.  Furthermore,  the  recirculation  zone  contains  three 
vortices.  This  is  physically  plausible  since  the  ramp  elongates  the  recirculation 
and  causes  a  small  vortex  to  be  snipped  off,  and  kinematic  considerations  require 
there  to  be  an  odd  number  of  vortices. 

Of  major  interest  is  a  comparison  of  the  velocity  profiles  in  the  shear 
layer  and  on  the  ramp.  These  results  are  shown  in  Figs.  28  and  29.  There  is 
fairly  good  agreement  between  the  computations  and  the  experiment.  In  particular 
the  presence  of  the  recompression  shock  on  the  ramp  is  evident  in  profile  F. 
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The  discrepancies  that  are  noted  are  due  in  part  to  the  mixing  length  turbulence 
model  that  was  employed.  Grid  refinement  on  the  ramp  and  moving  the  exit  plane 
further  downstream  could  also  aid  in  obtaining  better  agreement. 
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CONCLUSIONS 


In  this  report  the  development  of  a  Navier-Stokes  solution  procedure 
that  could  be  used  to  compute  the  aft  end  flow  field  over  missile  type  bodies 
is  presented.  Special  attention  was  given  to  handling  the  L-shaped  domain  that 
is  encountered,  the  treatment  of  the  sharp  reentrant  corner  and  the  application 
of  the  appropriate  boundary  conditions. 

An  adaptive  grid  option  was  exercised  for  a  model  back  step  problem 
where  the  free  shear  layer  was  tracked  and  the  grid  distribution  was 
automatically  adjusted  about  this  mean  free  shear  layer  curve.  Calculations 
were  performed  for  the  supersonic  turbulent  flow  over  a  10°  nozzle  boattail 
configuration  with  and  without  jet  exhaust.  The  qualitative  physical  behavior 
was  obtained  both  on  the  shoulder  and  in  the  wake.  Discrepancies  in  the  shoulder 
0^  signature  were  traced  at  least  in  part  to  the  turbulence  modeling.  For  the 
jet  off  case,  the  assumption  of  a  solid  base  rather  than  an  "open  nozzle  exit" 
could  account  for  some  of  the  observed  differences.  For  the  jet  on  case  greater 
grid  resolution  is  required  in  the  vicinity  of  the  corner. 

The  calculation  of  the  supersonic  turbulent  flow  over  a  2-D  back  step 
with  shear  layer  reattachment  on  a  20°  inclined  ramp  was  also  presented  and 
the  results  were  compared  with  experimental  data.  The  flow  field  results 
showed  qualitative  physical  behavior  while  the  velocity  profiles  in  the  shear 
layer  and  on  the  inclined  ramp  compared  well  with  data.  The  discrepancies 
that  were  observed  can  be  traced  again  in  part  to  turbulence  modeling,  grid 
resolution  and  placement  of  the  downstream  boundary. 

Further  evaluation  of  an  appropriate  turbulence  model  for  the  cases 
considered  is  crucial  in  obtaining  more  accurate  results.  Additional  investi¬ 
gation  of  the  effect  of  grid  distribution  and  resolution  is  also  warranted. 
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Figure  1  -  Schematic  of  Plume  Induced  Separation. 


F 

(a)  Before  Grid  Movement 
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(b)  After  Grid  Movement. 


Figure  3  -  Schematic  of  Adaptive  Grid 


Vector  Plot 


Figure  5  -  Comparison  of  Streamwise  Velocity  Profiles 
During  Forced  Grid  Movement. 
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Figure  6  -  Comparison  of  Fixed  Mesh  and  Adaptive  Mesh  Calculations, 
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Stations  and  Dimensions  in  Inches 
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4.075 

4.036 

3.993 

3.949 

3.905 

3.858 

3.809 

11.322 

11.569 

11.815 

12.062 

12 . 308 
12.555 
12.801 
13.048 
13.294 
13.541 
13.787 
14.034 
14.280 
14.527 
14.773 
15.020 
15.266 
15.513 
15.636 
15.759 
16.006 
16.252 
16.499 

3.760 

3.707 

3.653 

3.592 

3.530 

3.464 

3.397 

3.326 

3.252 

3.173 

3.089 

3.002 
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2.214 

2.125 

2.046 

Figure  7  -  Nozzle  Boattail  Geometry  (from  Ref.  5). 
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Coordinate 


A  Bottom  Surface  ^ 

D  Top  Surface  [  Data  of  Gausher 


Pressure  Distribution  on  Shoulder  of  Boattail  -  "Upstream  Boundary/Layer 


Vicinity  of  Base. 


a  Bottom  Surface  ^ 

□  Top  Surface  [  Data  of  Galigher  (Ref.  5) 


Pressure  Distribution  on  Shoulder  of  Boattail  -  Jet  Off 


Figure  17  -  Effect  of  Length  Scale  in  Wake  Region  on  Boattail 
Corner  Pressure  Jet  Off  Case. 
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Figure  22  -  Pressure  Distribution  on  Shoulder  of  Boattail  -  Jet  On. 


Figure  24  -  Schematic  or  Coordinate  Construction 


Computed  Contour  Plot  of  Pressure. 


•0381m  I  .0635m  I  .0889m 


SHEAR  LAYER  STATION 


0686m 


Figure  29  -  Velocity  Profiles  on  Ramp. 
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APPENDIX  A 


The  governing  conservation  equations  in  cylindrical-polar  coordinates  are 
transformed  using  the  Jacobian  transformation. 


yj  =  yj  (x, ,  x2,  x3,t) 
t  =  t 


(A-l) 


where  (x^,  x^)  =  (r,  0,  2).  The  resulting  equations  may  be  written  in  the 

following  compact  form: 


d(  JW) 


dr 


=  ~J{  a7<Jyi''*1  +  i|f/3iayrlJy'|F,) 


“£|  4l(jyi,iGi)]}  +  JS  +  JC 


where 


j  _  dy* 
y>  =  *T 


i.  _  dy 1 


y,i  = 


dx: 


Further,  the  coefficients  Bi>  y±,  are  given  by 


$1  =  r  *  z  r  1  ^3  =  * 


=  1  •  r2  =  —  •  v 3  =  1 


C  =  —  C 

=»i  rm  »  bj 


>  =1 


(A-2) 


(A-3) 


(A-4) 


and  m  ■  1  for  all  equations  except  the  x_-direction  momentum  equations  for  which 


The  vector  variables  used  in  Eq.  (A-2)  are  defined  as 


pUt 

pU.Uj 

PU  2 

/>u  2U, 

P  U3 

7;  =  rn 

p^\ 

P 

P  U  i 

pb 

P  hUj 

_Pk  _ 

- 1 

3~ 

(A-5)  (A-6) 


where  n  *  1  for  i  =  1  and  n  =  0  for  i  »  2,  3. 
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rTll 

Til 

r2T.2 

T \Z 

r  tI3 

Ti3 

0 

G,  = 

0 

-  rq. 

-  Qi 

iii  Y  t 

^  r, k.^ 
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m 

JC 

X 

alb" 

_ 1 

CA-7) 

(A-8) 

for  i  =2,3 


Note  that  the  velocity  components  (U^,  Uj,  U^)  are  the  cylindrical-polar  velocity 
components  written  in  cylindrical-polar  coordinates.  The  molecular  and  turbulent 
stress  tensors,  Eqs.  (6)  and  (8),  may  be  written  as 

-  i'1.f>(7-'iTl8ij  +  T<KBv-‘T-/,k>sii  <A-9) 


and  the  rate  of  strain  tensor  components  in  cylindrical-polar  coordinates  are 


Dm  = 


aui 

<3x, 


CL  -  ±  4^2.  +  ^L 

■22  r  dx2  r 


D,,= 


au3 

dx3 

_L  r 

r  d  ( 

U2  1 

2  L 

dx,  V 

r  ‘ 

j_  r 

au3 

aui 

2  L 

ax,  + 

d  x2 

_l  r 

j_  aua 

+ 

2  L 

r  dx2 

(A-10) 


and 


{rU, ) 


+ 


7 


dl>2 

dx2 


+ 


au3 
ax  3 


(A-ll) 


The  derivatives  required  in  Eqs.  (A-10  and  A-ll)  must  be  expressed  in  terms  of  the 
computational  coordinates  using  the  chain  rule,  Eq.  (34). 

Finally,  the  vector  £  contains  source  terms  and  certain  differential  terms 
which  do  not  conform  to  the  basic  structure  of  Eq.  (A-2) ,  and  the  vector  £  contains 
the  additional  curvature  terms  due  to  the  cylindrical-polar  coordinate  system. 


S  = 


0 

0 

0 

0 


(A-12) 
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APPENDIX  B 


Linearization  Technique 

A  number  of  techniques  have  been  used  for  implicit  solution  of  the  following 
first-order  nonlinear  scalar  equation  in  one  dependent  variable  4>(x,t): 


=F{<£)  <3G(£)/dx 


(Bl) 


Special  cases  of  Eq.  (Bl)  include  the  conservation  form  if  F(<J>)  =  1,  and  quasi- 
linear  flow  if  G($)  =  <)>.  Previous  implicit  methods  for  Eq.  (Bl)  which  employ 
nonlinear  difference  equations  and  also  methods  based  on  two-step  predictor- 
corrector  schemes  are  discussed  by  Ames  (Ref.  64,  p.  82)  and  von  Rosenburg 
(Ref.  65),  p.  56).  One  such  method  is  to  difference  nonlinear  terms  directly 
at  the  implicit  time  level  to  obtain  nonlinear  implicit  difference  equations; 
these  are  then  solved  iteratively  by  a  procedure  such  as  Newton's  method. 
Although  otherwise  attractive,  there  may  be  difficulty  with  convergence  in  the 
iterative  solution  of  the  nonlinear  difference  equations,  and  some  efficiency 
is  sacrificed  by  the  need  for  iteration.  An  implicit  predictor-corrector 
technique  has  been  devised  by  Douglas  and  Jones  (Ref.  66)  which  is  applicable 
to  the  quasilinear  case  (G  =  4> )  of  Eq.  (Bl) .  The  first  step  of  their  procedure 
is  to  linearize  the  equation  by  evaluating  the  non-linear  coefficient  as  F(<J>  )  - 
and  to  predict  values  of  using  either  the  backward  difference  or  the 

Crank-Nicolson  scheme.  Values  for  are  then  computed  in  a  similar  manner 

using  F(<J>n+^2)  and  the  Crank-Nicolson  scheme.  Gourlay  and  Morris  (Ref.  67) 
have  also  proposed  implicit  predictor-corrector  techniques  which  can  be  applied 
to  Eq.  (Al) .  In  the  conservative  case  (F=l) ,  their  technique  is  to  define 
G($)  by  the  relation  G(<l>)  ■  <j>G($)  when  such  a  definition  exists,  and  to  evaluate 
G(4>n+^)  using  the  values  for  4>n+  computed  by  an  explicit  predictor  scheme. 

With  G  thereby  known  at  the  implicit  time  level,  the  equation  can  be  treated  as 
linear  and  corrected  values  of  are  computed  by  the  Crank-Nicolson  scheme. 

A  technique  is  described  here  for  deriving  linear  implicit  difference 
approximations  for  nonlinear  differential  equations.  The  technique  is  based 
on  an  expansion  of  nonlinear  Implicit  terms  about  the  solution  at  the  known  time 
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level,  tn,  and  leads  to  a  one- step,  two-level  scheme  which,  being  linear  in 
unknown  (implicit)  quantities,  can  be  solved  efficiently  without  iteration. 
This  idea  was  applied  by  Richtmyer  and  Morton  (Ref.  68,  p.  203)  to  a  scalar 
nonlinear  diffusion  equation.  Here,  the  technique  is  developed  for  problems 
governed  by  £  nonlinear  equations  in  l  dependent  variables  which  are 
functions  of  time  and  space  coordinates.  Although  the  present  effort  concen¬ 
trates  upon  two  spatial  dimensions  and  time,  the  technique  will  be  described 
for  the  three-dimensional,  unsteady  equations. 


The  solution  domain  is  discretized  by  grid  points  having  equal  spacings 

12  3  12  3 

in  the  computational  coordinates,  Ay  ,  Ay  and  Ay  in  the  y  ,  y  and  y 

directions,  respectively,  and  an  arbitrary  time  step.  At.  The  subscripts 

12  3 

i,  j,  k  and  superscript  n  are  grid  point  indices  associated  with  y  ,  y  ,  y 

and  t,  respectively,  and  thus  4^  .  ,  denotes  ^(y'f,  y^,  y^,  tn) .  It  is 

l ,  J ,  k  i  k 

assumed  that  the  solution  is  known  at  the  n  level,  t  ,  and  is  desired  at 
the  (n+1)  level,  tn+^.  At  the  risk  of  an  occasional  ambiguity,  one  or  more 
of  the  subscripts  is  frequently  omitted,  so  that  4°  is  equivalent  to  .  ,  . 

i,  J  ,K 

Although  present  attention  is  focused  on  the  compressible  Navier-Stokes 
equations,  the  numerical  method  employed  is  quite  general  and  is  formally 
derived  for  systems  of  governing  equations  which  have  the  following  form: 


dH[<p)/d\  =  2>(<£)  +S(<£)  (B2) 

where  <t>  is  a  column  vector  containing  £  dependent  variables,  H  and  S  are 
column  vector  functions  of  $,  and  2>  is  a  column  vector  whose  elements  are 
spatial  differential  operators  which  may  be  multidimensional.  The  generality 
of  Eq.  (B2)  allows  the  method  to  be  developed  concisely  and  permits  various 
extensions  and  modifications  (e.g.,  noncartesian  coordinate  systems,  turbulence 
models)  to  be  made  more  or  less  routinely.  It  should  be  emphasized,  however, 
that  the  Jacobian  3H/'d«j>  must  usually  be  nonsingular  if  the  ADI  techniques  as 
applied  to  Eq.  (B2)  are  to  be  valid.  A  necessary  condition  is  that  each 
dependent  variable  appear  in  one  or  more  of  the  governing  equations  as  a  time 
d*  i^ative.  An  exception  would  occur  if  for  instance,  a  variable  having  no  time 
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derivative  also  appeared  in  only  one  equation,  so  that  this  equation  could  be 
decoupled  from  the  remaining  equations  and  solved  a  posteriori  by  an  alternate 
method.  As  a  consequence,  the  present  method  is  not  directly  applicable  to  the 
incompressible  Navier-Stokes  equations  except  in  one-dimension,  where  ADI 
techniques  are  unnecessary.  For  example,  the  velocity-pressure  form  of  the 
incompressible  equations  has  no  time  derivative  of  pressure,  whereas  the 
vorticity-stream-function  form  has  no  time  derivative  of  stream  function.  For 
computing  steady  solutions,  however,  the  addition  of  suitable  "artificial"  time 
derivatives  to  the  incompressible  equations,  as  was  done  in  Chorin's  (Ref.  69) 
artificial  compressibility  method,  would  permit  the  application  of  the  present 
method.  Alternatively,  a  low  Mach  number  solution  of  the  compressible  equations 
can  be  computed. 

The  linearized  difference  approximation  is  derived  from  the  following 
implicit  time-difference  replacement  of  Eq.  (B2) : 

(Hn+’-Hn)/At=)9[2>(^n+,)+Sn  +  ,]  +  (l-/3)[^(<#)n)+Sn]  (B3) 


where,  for  example,  Hn+1  =  H((f>n+1 ) .  The  form  of  -2>  and  the  spatial  differencing 
are  as  yet  unspecified.  A  parameter  0(0  <  6  <  1)  has  been  introduced  so  as  to 
permit  a  variable  centering  of  the  scheme  in  time.  Equation  (B3)  produces  a 
backward  difference  formulation  for  6=1  and  a  Crank-Nicolson  formulation  for 
B  =  1/2. 

The  linearization  is  perfcrmed  by  a  two-step  process  of  expansion  about  the 
known  time  level  tn  and  subsequent  approximation  of  the  quantity  (3<|>/3t)nAt, 
which  arises  from  chain  rule  differentiation,  by  (<j>n+^-<)>n)  .  The  result  is 

Hn+I  =  Hn  +  UH/<J(£)n  (<£n+‘-<£n)  +  0(Af)2  (B4a) 

Sn+l=Sn+(dS/<fc£)n  (4>0  +  '-4>n)  +0(  At)2  (B4b) 
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(B4c) 


2>(<£n+')  =  -2>(^>n)+(a  ^>/a<^)(^n+,-0n)  +  o(At)5 


The  matrices  3H/3<}>  and  3S/34>  are  standard  Jacobians  whose  elements  are 

defined,  for  example,  by  (3H/3<j>)  =  3H  /3<f>  .  The  operator  elements  of  the 

qr  q  r 

matrix  3.2>/3*  are  similarly  ordered,  i.e.,  (3.2>/3<j>)  H3.2>/3$  ;  however,  the 

cjr  <1  r 

intended  meaning  of  the  operator  elements  requires  some  clarification.  For  the 
qtl*  row,  the  operation  /3$)n  (4>n+^  -  <t>n)  is  understood  to  mean  that 

q  n 

(3/3t-£  [$(x,y,z,t)  ]  }n  At  is  computed  and  that  all  occurrences  of  ( 3<#>r/3t) 
arising  from  chain  rule  differentiation  are  replaced  by  -  <I>”)/At. 

After  linearization  as  in  Eqs.  (B4) ,  Eq.  (A3)  becomes  the  following  linear 
implicit  time-differenced  scheme: 


(dHn/dWn  +  l 1  —  «#> n) /At  =2>(<£n)  +Sn  +  (3  (d2>/d<f>  +  dSn/d<f>H<t>n+'-<f>n)  (B5) 


Although  Hn+1  is  linearized  to  second  order  in  Eq.  (B4)  ,  the  division  by  At  in 
Eq.  (B3)  introduces  an  error  term  of  order  At.  A  technique  for  maintaining 
formal  second-order  accuracy  in  the  presence  of  nonlinear  time  derivatives  is 
discussed  by  McDonald  and  Briley  (Ref.  32) ,  however,  a  three-level  scheme 
results.  Second-order  temporal  accuracy  can  also  be  obtained  (for  0  =  1/2)  by 
a  change  in  dependent  variable  to  $  =  H (<p) ,  provided  this  is  convenient,  since 
the  nonlinear  time  derivative  is  then  eliminated.  The  temporal  accuracy  is 
independent  of  the  spatial  accuracy. 

On  examination,  it  can  be  seen  that  Eq.  (B5)  is  linear  in  the  quantity 
($n+^  -  $n)  and  that  all  other  quantities  are  either  known  or  evaluated  at  the 
n  level.  Computationally,  it  is  convenient  to  solve  Eq.  (B5)  for  (<f>n+^  -  $n) 
rather  than  <J>  .  This  both  simplifies  Eq.  (B5)  and  reduces  roundoff  errors, 

since  it  is  presumably  better  to  compute  a  small  O(At)  change  in  an  0(1) 
quantity  than  the  quantity  itself.  To  simplify  the  notation,  a  new  dependent 
variable  \jj  defined  by 

*  =  4>-4>n  (B6) 
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1 


<  ■-..■•■-'■w+H  l  m  >i^ 


is  introduced,  and  thus  =  4>  n+*  -  41°,  and  t|?n  =  0.  It  is  also  convenient 

to  rewrite  Eq.  (A5)  in  the  following  simplified  form: 

(A+At^)^n  +  I  =  At  [$  (^n)+Sn]  (B7a) 

where  the  following  symbols  have  been  introduced  to  simplify  the  notation: 


A  =  dHn/d4>  -j3At(dSn/d<£) 


(B7b) 


1=  -p{d  2>/d<f>) 


(B7c) 


It  is  noted  that  is  a  linear  transformation  and  thus  ./(0)  *  0.  Further¬ 

more,  if  ,2>(<j>)  is  linear,  then  I  (41)  =  -B  .2>(if))  - 

Spatial  differencing  of  Eq.  (B7a)  is  accomplished  simply  by  replacing 

derivative  operators  such  as  3/3yi,  3^/3yi3yi  by  corresponding  finite  differ- 

2 

ence  operators,  ,  D^.  Henceforth,  it  is  assumed  that.  and  I  have  been 

discretized  in  this  manner,  unless  otherwise  noted. 

Before  proceeding,  some  general  observations  seem  appropriate.  The  fore¬ 
going  linearization  technique  assumes  only  Taylor  expandability,  an  assumption 
already  implicit  in  the  use  of  a  finite  difference  method.  The  governing 
equations  and  boundary  conditions  are  addressed  directly  as  a  system  of  coupled 
nonlinear  equations  which  collectively  determine  the  solution.  The  approach 
thus  seems  more  natural  than  that  of  making  ad  hoc  linearization  and  decoupling 
approximations,  as  is  often  done  in  applying  implicit  schemes  to  coupled  and/or 
nonlinear  partial  differential  equations.  With  the  present  approach,  it  is  not 
necessary  to  associate  each  governing  equation  and  boundary  condition  with  a 
particular  dependent  variable  and  then  to  identify  various  "nonlinear 
coefficients"  and  "coupling  terms"  which  must  then  be  treated  by  lagging, 
predictor-corrector  techniques,  or  iteration.  The  Taylor  expansion  procedure 
is  analogous  to  that  used  in  the  generalized  Newton-Raphson  or  quasi-linearization 
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methods  for  iterative  solution  of  nonlinear  systems  by  expansion  about  a  known 
current  guess  at  the  solution  (e.g.,  Bellman  &  Kalaba,  Ref.  70).  However,  the 
concept  of  expanding  about  the  previous  time  level  apparently  had  not  been 
employed  to  produce  a  noniterative  implicit  time-dependent  scheme  for  coupled 
equations,  wherein  nonlinear  terms  are  approximated  to  a  level  of  accuracy 
commensurate  with  that  of  the  time  differencing.  The  linearization  technique 
also  permits  the  implicit  treatment  of  coupled  nonlinear  boundary  conditions, 
such  as  stagnation  pressure  and  enthalpy  at  subsonic  inlet  boundaries,  and  in 
practice,  this  latter  feature  was  found  to  be  crucial  to  the  stability  of  the 
overall  method  (Ref .  32 ) . 

Application  of  Alternating-Direction  Techniques 

Solution  of  Eq.  (B7a)  is  accomplished  by  application  of  an  alternating- 
direction  implicit  (ADI)  technique  for  parabolic-hyperbolic  equations.  The 
original  ADI  method  was  introduced  by  Peaceman  and  Rachford  (Ref.  71)  and 
Douglas  (Ref.  72);  however,  the  alternating-direction  concept  has  since  been 
expanded  and  generalized.  A  discussion  of  various  alternating-direction 
techniques  is  given  by  Mitchell  (Ref.  73)  and  Yanenko  (Ref.  74  ). 

The  present  technique  is  simply  an  application  of  the  very  general 
procedure  developed  by  Douglas  and  Gunn  (Ref.  35)  for  generating  ADI  schemes 
as  perturbations  of  fundamental  implicit  difference  schemes  such  as  the 
backward-difference  or  Crank-Nicolson  schemes. 

For  the  present,  it  will  be  assumed  that  (<f> )  contains  derivatives  of 

12  3 

first  and  second  order  with  respect  to  y  ,  y  and  y  ,  but  no  mixed  derivatives. 

In  this  case,  can  be  split  into  three  operators,  associated 

12  3  l  Z  3 

with  the  y  ,  y  and  y  coordinates  and  each  having  the  functional  form 

.2^  *  S/Sy*,  3^/3y*3y*).  Equation  (A7a)  then  becomes 

[a  +  &!(7,  +  +./3)]*n'H=  At[(2>,  +^2  +  ^3)^n  +  Sn]  (B8) 

Recalling  that  ,/(^n)  =  0,  the  Douglas-Gunn  representation  of  Eq.  (B8)  can  be 
written  as  the  following  three-step  solution  procedure: 
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(A+At  Jx  )f*=  At[(  2>,  +^2  +  2>3)*n+Sn] 


(B9a) 


(A  +  At  J2W**  =Af* 


(B9b) 


(A  +  At  Jz)^  +  l  =  A+** 


(B9c) 


where  iJi  and  ip  are  intermediate  solutions.  It  will  be  shown  subsequently 

that  each  of  Eqs.  (B9)  can  be  written  in  narrow  block-banded  matrix  form  and 

*  ** 

solved  by  efficient  block-elimination  methods.  If  i/;  and  if*  are  eliminated, 
Eqs.  (B9)  become 


(A  +  At4)A"'(A  +  AJ  J2  ) A- 1  ( A  +  At  )V'"'H  =  At 


[cV2>2  +  3>  3)  <#>n+sn] <B10> 


If  the  multiplication  on  the  left-hand  side  of  Eq.  (BIO)  is  performed,  it  becomes 

2 

apparent  that  Eq.  (BIO)  approximates  Eq.  (B8)  to  order  (At)  .  Although  the 
stability  of  Eqs.  (B9)  has  not  been  established  in  circumstances  sufficiently 
general  to  encompass  the  Navier-Stokes  equations,  it  is  often  suggested  (e.g., 
Richtmyer  &  Morton,  Ref.  68,  p.  215)  that  the  scheme  is  stable  and  accurate  under 
conditions  more  general  than  those  for  which  rigorous  proofs  are  available. 

This  latter  notion  was  adopted  here  as  a  working  hypothesis  supported  by  favor¬ 
able  results  obtained  in  actual  computations  (e.g..  Refs.  33  and  34). 

A  major  attraction  of  the  Douglas-Gunn  scheme  is  that  the  intermediate 

*  **  n+l 

solutions  and  <J<  are  consistent  approximations  to  <l>  .  Furthermore, 

,n  *  **  n+l 

for  steady  solutions,  ^  ^  ^  independent  of  At.  Thus,  physical 

n+l 

boundary  conditions  for  \f>  can  be  used  in  the  intermediate  steps  without  a 
serious  loss  in  accuracy  and  with  no  loss  for  steady  solutions.  In  this  respect, 
the  Douglas-Gunn  scheme  appears  to  have  an  advantage  over  locally  one-dimensional 
(LOD)  or  "splitting"  schemes,  and  other  schemes  whose  intermediate  steps  do  not  fl 

satisfy  the  consistency  condition.  The  lack  of  consistency  in  the  intermediate  1 

steps  complicates  the  treatment  of  boundary  conditions  and,  according  to  Yanenko  1 


(Ref.  74,  p.  33),  does  not  permit  the  use  of  asymptotically  large  time  steps. 

It  is  not  clear  that  this  advantage  of  the  Douglas-Gunn  scheme  would  always 
outweigh  other  benefits  which  might  be  derived  from  an  alternative  scheme. 
However,  since  the  ADI  scheme  can  be  viewed  as  an  approximate  technique  for 
solving  the  fundamental  difference  scheme,  Eq.  (B7a) ,  alternate  techniques  can 
readily  be  used  within  the  present  formulation. 

It  is  worth  noting  that  the  operator  lb  can  be  split  into  any  number  of 
components  which  need  not  be  associated  with  a  particular  coordinate  direction. 
As  pointed  out  by  Douglas  and  Gunn  (Ref.  41),  the  criterion  for  identifying 
sub-operators  is  that  the  associated  matrices  be  "easily  solved"  (i.e.,  narrow-' 
banded).  Thus,  mixed  derivatives  can  be  treated  implicitly  within  the  ADI 
framework,  although  this  would  increase  the  number  of  intermediate  steps  and 
thereby  complicate  the  solution  procedure.  Finally,  only  minor  changes  are 
introduced  if,  in  the  foregoing  development  of  the  numerical  method,  H,  .2),  and 
S  are  functions  of  the  spatial  coordinates  and  time,  as  well  as  <J> . 
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